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has  been  previously  studied  and  large  suppression  of  the  weak  signal  has  hppn— 
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obtained.  Even  though  preliminary  analytical  work  indicated  that  little  suppres' 
sion  is  possible  in  Case  II,  the  present  research  was  undertaken  to  determine 
the  amount  of  suppression  possible.  Various  exact  and  approximate  preshock 
solutions  were  reviewed  in  the  time  domain  (Earnshaw  solution)  and  frequency 
domain  (Fenlon  solution).  Postshock  solutions  were  based  on  weak  shock  theory 
which  was  implemented  in  a computer  algorithm,  corrections  for  attenuation  and 
dispersion  were  incorporated.  Predictions  thus  obtained  were  compared  to  expe- 
imental  results.  The  experiments,  performed  in  a 30  m progressive  wave  tube, 
were  restricted  to  Case  II.. — For  the  weak  signal,  the  range  of  source  SPL  was 
104  to  121  dB  and  the  range  ofvfrequency  was  0.6  to  1 kHz.  The  corresponding 
ranges  for  the  pump  were  125  to  OB~riB--and^l . 5 to  6.6  kHz,  respectively.  The 
measurement  distances  ranged  from  3.9  m to  TTrn7'^The  agreement  between  exper- 
imental results  and  theoretical  predictions  was  excellent.  In  particular  the 
experiments  confirmed  the  theoretical  prediction  that  little  suppression  of  the 
weak  signal  is  obtained  for  Case  II  either  before  or  after  shock  formation. 

The  sidebands  in  both  cases  form  around  the  high  frequency  and  its  harmonics. 

Thus  another  general  conclusion  is  that  the  low  frequency  signal  modulates  the 
high  frequency  signal,  regardless  of  their  respective  amplitudes. 


UNCLASSIFIED 

SECuRITYTCLASSIFldATlON  OF  THIS  PAOEfWhin  Dele  Entered) 


ARL-TR-77-22 

May  1977 


THE  SUPPRESSION  OF  SOUND  BY  SOUND  OF  HIGHER  FREQUENCY 


William  L.  Willshire,  Jr. 


Office  of  Naval  Research 
Contract  N00014 - 75 - C - 0867^ 

Air  Force  Office  of  Scientific  Research 
Contract  F44620  - 76  - C - 0040 


APPLIED  RESEARCH  LABORATORIES 

THE  UNIVERSITY  OF  TEXAS  AT  AUSTIN 

AUSTIN,  TEXAS  7871 2 


APPROVED  FOR  PUBLIC 
RELEASE.  DISTRIBUTION 
UNLIMITED 


FOREWORD 


| 


This  report  is  an  adaptation  of  William  L.  Willshire's  thesis, 
"The  Suppression  of  Sound  by  Sound  of  Higher  Frequency,"  which  was 
written  for  the  degree  of  Master  of  Science  in  Engineering  at  The 
University  of  Texas  at  Austin.  Mr.  Willshire  was  enrolled  in  the 
Department  of  Mechanical  Engineering,  and  the  date  of  his  degree  is 
May  1977- 

The  research  was  carried  out  at  the  Applied  Research  Laboratories 
and  was  supported  by  the  Office  of  Naval  Research  under  Contract 
N00014-75-C-0867  and  by  the  Air  Force  Office  of  Scientific  Research 
under  Contract  F44620-76-C-0040.  Technical  monitors  were 
Dr.  Walter  Madigosky  and  Dr.  Logan  E.  Hargrove  for  ONR,  and 
Lt.Col.  R.  C.  Smith  and  Lt.Col.  L.  W.  Ormand  for  APOSR. 

David  T.  Blackstock 
Supervisor 


ABSTRACT 


The  purpose  of  this  research  was  to  study  the  suppression  of  one 
sound  by  another  sound  of  higher  frequency.  More  specifically,  the  planar 
propagation  in  air  of  a finite  amplitude  tone  (the  pump)  and  a smaller 
amplitude  tone  (the  weak  signal)  of  lower  frequency  was  investigated  both 
theoretically  and  experimentally.  Suppression  is  caused  by  a modulation 
of  the  weak  signal  by  the  pump.  The  acoustic  energy  of  the  weak  signal  is 
not  lost  (absorbed)  but  rather  redistributed  to  other  parts  of  the  spec- 
trum. The  modulation  process  leads  to  the  formation  of  sidebands  about 
the  high  frequency  signal.  The  amount  of  suppression  depends  qualita- 
tively on  whether  the  pump  frequency  is  lower  than  the  weak  signal 
frequency  (Case  I)  or  vice  versa  (Case  II).  Case  I has  been  previously 
studied  and  large  suppression  of  the  weak  signal  has  been  obtained.  Even 
though  preliminary  analytical  work  indicated  that  little  suppression  is 
possible  in  Case  II,  this  research  was  undertaken  to  determine  the  amount 
of  suppression  possible. 

The  theoretical  investigation  was  concerned  with  both  Case  I 
and  Case  II;  a general  notation  was  adopted  so  that  the  analysis  of  a 
single  solution  was  applicable  to  both  cases.  Various  exact  and  approxi- 
mate preshock  solutions  were  reviewed  in  the  time  (Earnshaw  solution)  and 
frequency  (Fenlon  solution)  domains.  These  solutions  were  then  compared 
with  each  other  and  ranked  with  regard  to  their  exactness.  Because  general 
solutions  valid  in  the  postshock  region  were  not  available,  postshock 
solutions  were  developed  based  on  weak  shock  theory.  Weak  shock  theory  was 
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implemented  in  a computer  algorithm  with  corrections  for  attenuation  and 
dispersion.  Predictions  obtained  from  the  use  of  the  program  were  com- 
pared to  experimental  results.  The  experimental  work  was  restricted  to 
Case  II.  For  the  weak  signal  the  experiments  (performed  in  a 30  m 
progressive  wave  tube)  covered  a range  of  source  SPL  of  10 4 to  121  dB  and 
a frequency  range  of  0.6  to  1 kHz.  The  corresponding  ranges  for  the  pump 
were  125  to  158  dB  and  1.5  to  6.6  kHz,  respectively.  The  measurement 
distances  ranged  from  3*9  m to  27  m.  The  agreement  between  experimental 
results  and  theoretical  predictions  was  excellent.  In  particular  the 
experiments  confirmed  the  theoretical  prediction  that  little  suppression 
of  the  weak  signal  is  obtained  for  Case  II  either  before  or  after  shock 
formation.  The  sidebands  in  both  cases  form  around  the  high  frequency 
signal  and  the  energy  for  the  sideband  formation  comes  primarily  from  the 
high  frequency  signal  (including  its  harmonics).  Thus  another  general 
conclusion  is  that  the  low  frequency  signal  modulates  the  high  frequency 
signal,  regardless  of  their  respective  amplitudes. 
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INTRODUCTION  AND  LITERATURE  REVIEW 

The  nonlinear  interaction  of  two  or  more  acoustic  signals  has 
been  a very  active  area  of  study  in  the  past  decade  or  so.  In  a note- 
worthy early  investigation  Thuras,  Jenkins,  and  O'Neil'*'  (1935)  showed  that 

when  two  intense  tones  propagate  in  a plane  wave  tube,  sum  and  difference 

2 

frequency  components  are  created.  In  1956  Ingard  and  Pridmore -Brown 
reported  a study  of  the  interaction  of  two  intersecting  beams  of  intense 
sound  and  used  the  phrase  "the  scattering  of  sound  by  sound"  to  describe 
the  interaction.  In  a theoretical  investigation  Westervelt^  (1957) 
claimed  that  the  interaction  would  take  place  only  for  collinear  propaga- 
tion of  the  two  beams.  The  work  on  the  scattering  of  sound  by  sound  led 
tc  the  invention  of  the  parametric  array  by  Westervelt^  (i960).  In  turn, 
a recent  spinoff  of  research  on  the  parametric  array  has  been  the 

c 

"absorption  of  sound  by  sound,"  a term  coined  by  Westervelt^  (1973)  to 
describe  the  use  of  one  sound  to  decrease  the  amplitude  of  a second  sound. 

A more  descriptive  phrase  for  the  phenomenon  is  "suppression  of  sound  by 
sound"  because  the  energy  lost  by  the  second  sound  is  not  absorbed  but 
rather  redistributed  to  other  frequencies.  For  an  excellent  review  of 
works  related,  both  directly  and  indirectly,  to  suppression  of  sound  by 
sound  see  Schaffer^  (1975) • Suppression  of  sound  by  sound  is  the  subject 
of  this  thesis.  More  specifically,  the  subject  of  this  thesis  is  the 
suppression  of  a small  signal  by  a higher  frequency  intense  signal. 

1 
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An  underwater  experiment  with  spherical  waves  performed  by 
Moffett,  Konrad,  and  Corcella  (1972)  demonstrated  the  phenomenon  of 
suppression  of  sound  by  sound.  The  experiment  was  intended  to  test  a 
parametric  array  whose  primary  frequencies  were  widely  separated. 

Moffett  et  al.  found  that  the  amplitude  of  the  high  frequency  primary 
was  decreased  by  17  dB  when  the  low  frequency  primary  was  turned  on.  On 
the  other  hand,  little  decrease  in  the  low  frequency  primary  was  observed 
when  the  high  frequency  primary  was  turned  on.  The  explanation  given  for 
the  decrease  in  the  first  case  was  that  the  high  frequency  primary  under- 
went a frequency  modulation  by  the  low  frequency  primary.  The  modulation 
caused  a production  of  sidebands  about  the  high  frequency  "carrier"  - at 
the  carrier's  expense.  When  the  roles  of  the  two  primaries  were  reversed 
(the  second  case),  a modulation  still  occurred,  but  one  that  had  little 
effect  on  the  amplitude  of  the  low  frequency  wave. 

The  experiments  of  Moffett  et  al.  show  that  two  qualitatively 

different  cases  are  to  be  considered  when  a finite  amplitude  tone  is  used 

* 

to  suppress  a smaller  amplitude  tone.  In  Ca6e  I the  high  intensity  wave, 
which  will  be  called  the  pump,  is  lower  in  frequency  than  the  weak  signal. 
In  Case  II  the  pump  is  higher  in  frequency  than  the  weak  signal. 

The  interaction  processes  in  Case  I and  Case  II  are  graphically 
illustrated  in  Fig.  1.1  and  Fig.  1.2,  respectively.  Here  u (vertical 
axis)  is  particle  velocity,  t is  time,  x is  propagation  distance,  cq  i6 
the  small  signal  propagation  velocity,  t'=t-(x/cQ)  is  retarded  time,  cu  is 
angular  frequency,  and  x is  shock  formation  distance.  In  general  the 

Moffett  et  al.  studied  the  interaction  of  two  intense  signals,  but  the 
distinction  between  the  two  cases  may  still  be  made. 
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subscript  p refers  to  the  pump  and  w to  the  weak  signal.  The  combined 
signal  is  represented  by  u,  u ' is  the  pump  propagated  independently, 
that  is,  in  the  absence  of  the  weak  signal,  and  uw'  is  the  weak  signal 
propagated  independently.  A zero  preceding  a subscript  refers  to  source 

g 

amplitude.  Lossless  theory  has  been  used  to  compute  the  waveforms. 

First  consider  Case  I.  The  low  frequency  pump  by  itself  is  shown  in 
Figs,  l.l(a)  and  (b),  first  at  the  source  (x=0)  and  then  after  some  finite 
amplitude  distortion  has  occurred  (x=Xp/2).  Figures  l.l(c)  and  (d)  show 
the  weak  signal  by  itself  at  the  same  propagation  distances.  The  combined 
signal  at  the  same  distances  is  given  in  Figs,  l.l(e)  and  (f).  The  dis- 
torted pump  waveform  [1.1(b)]  is  subtracted  from  the  distorted  combined 
waveform  [1.1(f)]  to  illustrate  the  effect  the  pump  distortion  has  on  the 
weak  high  frequency  signal.  The  result  is  shown  in  Fig.  l.l(g);  an 
amplification  factor  of  3 has  been  applied  to  emphasize  the  result.  Next 
consider  Case  II.  The  source  waveforms  for  the  weak  signal,  the  high 
frequency  pump,  and  the  combination  of  the  pump  and  weak  signal  are  shown 
in  Figs.  1.2(a),  (c),  and  (e),  respectively.  The  three  source  waveforms 
are  shown  propagated  separately  an  equal  distance  (x=x^)  in  Figs.  1.2(b), 
(d),  and  (f),  respectively.  Subtraction  of  the  distorted  waveform 
[1.2(d)1  from  the  distorted  combined  waveform  [1.2(f)]  leads  to  the  wave- 
form shown  in  Fig.  1.2(g).  Again,  an  amplification  factor  equal  to  3 has 
been  used  to  emphasize  the  result. 

Through  Figs.  1.1  and  1.2  a simple  qualitative  explanation  may 
be  given  for  the  two  different  observations  of  Moffett  et  al.  In  Case  I 
a frequency  modulation  of  the  weak  signal  is  clearly  visible  [1.1(g)] . 

The  frequency  modulation  is  successful  because  the  frequency  of  the  pump 
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is  less  than  that  of  the  weak  signal:  the  nonlinear  distortion  of  a 

single  pump  cycle  affects  many  cycles  of  the  weak  signal.  Little  varia- 
tion of  the  weak  signal  is  achieved  for  Case  II,  however.  A high 
frequency  ripple  at  approximately  the  frequency  of  the  pump  is  observed 
in  Fig.  1.2(g).  In  fact,  the  effect  most  apparent  in  the  Case  II  figure 
is  a low  frequency  modulation  of  the  pump  by  the  weak  signal  [Fig.  1.2(f)]. 
The  distortion  of  a single  pump  cycle  influences  only  a fraction  of  a weak 
signal  cycle.  The  ripple  on  the  weak  signal  is  the  result  of  the  dis- 
tortion of  successive  pump  cycles. 

Additional  insight  into  the  suppression  of  sound  by  sound 
interaction  is  obtained  by  viewing  the  process  in  the  frequency  domain. 
Figures  1.3  and  1.4  show  spectra  of  waveforms  similar  to  those  in 
Figs.  1.1  and  1.2.  The  spectra  are  actually  sample  data  taken  from 
experiments  that  are  described  in  Chapter  TV.  Figure  1.3  illustrates 
Case  I;  Fig.  1.4  illustrates  Case  II.  When  one  acoustic  signal  is 
suppressed  by  a second,  the  energy  taken  from  the  first  is  redistributed 
to  sidebands  centered  around  the  higher  frequency  wave,  and,  in  Case  II, 
its  harmonics.  In  Case  I the  pump  modulates  the  high  frequency  wave  and 
transfers  energy  from  the  weak  signal  to  sidebands  located  at  u^inco^ 

(n  an  integer).  Again  in  Case  II,  the  low  frequency  signal  modulates  the 
high  frequency  signal,  but  the  sideband  energy  (located  at  mu^±o^,  m an 
integer)  appears  to  be  drawn  from  the  pump  and  its  harmonics,  not  the  weak 
signal.  In  both  cases  the  energy  for  the  sideband  formation  primarily 
originates  from  the  high  frequency  signal.  In  Case  I the  energy  redis- 
tributed to  the  sidebands  from  the  weak  signal  substantially  decreases 
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the  amplitude  of  the  weak  signal;  in  Case  II  the  energy  redistributed  from 

the  pump  and  its  harmonics  has  little  effect  on  their  amplitudes. 

Although  qualitative  discussions  are  useful  to  achieve  a basic 

understanding  of  the  phenomenon  of  suppression  of  sound  by  sound,  further 

understanding  of  the  phenomenon  requires  a quantitative  investigation. 

A way  to  begin  the  investigation  is  to  use  a solution  by  Fenlon  for  dis- 

o 

tortion  of  sound  from  a two-frequency  source.  Although  his  solution  is 
limited  to  lossless  propagation  and  to  distances  less  than  x,  it  provides 
important  information  about  the  parameters  that  are  significant  and  about 
the  amount  of  suppression  that  may  be  obtained.  The  predicted  result  for 
the  amplitude  of  the  weak  signal  for  both  cases,  provided  the  source 
amplitude  of  the  weak  signal  is  much  less  than  that  of  the  pump,  is 


u , = u J (fta  ) 
weak  ow  o p 


(1.1) 


where  the  subscript  weak  refers  to  the  weak  signal  of  the  combined  wave- 
form, J is  the  zero  order  Bessel  function  of  the  first  kind,  is  the 
ratio  of  the  weak  signal  frequency  to  the  pump  frequency,  and  a^=x/x^. 
Complete  suppression  of  the  weak  signal  may  be  achieved  at  the  zeroes  of 
Jq,  the  first  one  occurring  at  fia^=2.4.  Here  a quantitative  difference 
between  Case  I and  Case  II  may  be  seen.  First  note  that  the  prediction  is 
valid  only  when  ct^<1.  Even  with  this  restriction,  however,  in  Case  I 
may  easily  be  chosen  large  enough  that  £2ct^=2.4.  Thus  a complete  null  may 
in  theory  be  achieved.  In  an  experimental  study  of  Case  I interaction, 
Schaffer^  observed  reductions  of  as  much  as  4o  dB  in  the  weak  signal  when 
he  varied  the  product  over  a range  including  the  null  value  2.4.  In 
Case  II,  on  the  other  hand,  fi  has  values  less  than  one.  Therefore,  the 
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largest  permissible  value  of  the  argument  of  the  Bessel  function  is  fia 


At  this  extreme  value,  the  reduction  of  the  weak  signal  is  only  2.4  dB. 


The  effectiveness  of  the  suppression  of  the  weak  signal  is  shown 

for  both  cases  in  Pig.  1.5-  Plot  (a)  is  for  Case  I;  the  horizontal  axis 

is  a and  the  vertical  axis  is  the  decrease  [20  log,„  J (Ha  )]  in  the  weak 
p B10  ox  p 

signal.  Three  different  values  of  fi  are  used.  When  fi=4,  the  range  of  the 
argument  of  the  Bessel  function  includes  its  first  zero,  and  total 
suppression  of  the  weak  signal  is  achieved.  Figure  1.5  plot  (b)  shows 
the  suppression  of  the  weak  signal  for  Case  II.  Three  different  values  of 
fi  are  used.  Within  the  limits  of  the  validity  of  the  theory  and  within 


possible  values  of  fi,  there  is  little  suppression  of  the  weak  signal. 


Other  quantitative  approaches  besides  Fenlon's  solution  have 
been  used  to  investigate  the  nonlinear  interaction  between  two  acoustic 
signals.  A frequency  modulation  ( PM) -type  solution  has  been  used  by 
Schaffer^  and  others  to  describe  the  Case  I interaction.  Schaffer  found 
that  PM  solution  and  Penlon  solution  give  the  same  prediction  for  the 
amplitude  of  the  weak  signal  but  not  for  the  amplitudes  of  the  sidebands. 
In  the  FM  solution  the  sidebands  are  symmetric  in  magnitude  about  the 
weak  signal.  In  the  Fenlon  solution,  however,  the  amplitudes  are 
asymmetrical  about  the  weak  signal;  the  upper  sidebands  have  the  larger 
amplitudes.  An  FM  solution  may  also  be  derived  for  the  Case  II  inter- 
action but,  as  with  the  Fenlon  solution,  the  FM  solution  is  valid  only 
for  shock  free  waveforms. 

Little  suppression  of  the  weak  signal  of  Case  II  takes  place 
before  shocks  form.  An  investigation  is  needed  to  determine  whether 
additional  suppression  of  the  weak  signal  can  be  achieved  after  shock 
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formation.  Little  analytical  work  has  been  done  on  the  interaction  of  two 
acoustic  signals  after  shock  formation.  Although  an  early  treatment  was 
given  by  Naugol'nykh,  Soluyan,  and  Khokhlov10  (1963),  and  Moffett  et  al.11 
(1977)  discuss  the  phenomenon,  no  general  solution  of  the  suppression 
problem  valid  after  shock  formation  has  been  found.  The  weak  shock  method 
is  an  available  mathematical  procedure  which  yields  valid  solutions  beyond 
the  formation  of  shocks.  A substantial  amount  of  the  theoretical  work 
involved  with  this  thesis  deals  with  the  application  of  weak  shock  theory 
to  investigate  Case  II  suppression  at  distances  x>x. 

Although  suppression  of  sound  by  sound  seems  to  be  much  more 
effective  in  Case  I than  in  Case  II,  at  least  for  the  shock-free  region, 
Case  II  is  generally  more  attractive  when  applications  are  considered. 

An  application  that  immediately  comes  to  mind  is  the  controlled  reduction 
of  acoustic  signals.  The  practicality  of  using  Case  I for  acoustic 
control  is  limited.  Although  much  reduction  of  a high  frequency  signal 
can  be  obtained,  the  frequency  region  about  the  reduced  signal  is  swamped 
by  harmonic  and  sideband  "pollution.”  Figure  1.3(c)  illustrates  the 
situation.  Considerable  reduction  of  the  weak  signal  has  been  achieved, 
but  much  acoustic  energy  is  located  at  nearby  frequencies.  The  applica- 
tion of  Case  I for  acoustic  reduction  is  limited  to  a narrow  frequency 
band,  in  particular,  a band  whose  width  is  no  greater  than  the  frequency 
of  the  pump.  Case  II,  on  the  other  hand,  has  more  possibilities; 

Fig.  l.^(c)  shows  why.  Because  the  pump  frequency  in  Case  II  is  greater 
than  that  of  the  weak  signal,  the  acoustic  pollution  is  higher  in  fre- 
quency than  the  weak  signal.  If  the  pump  can  reduce  the  lower  frequency 
weak  signal,  additional  acoustical  energy  will  be  required  that  is  not 
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close  in  frequency  to  the  weak  signal.  The  reduction  bandwidth  is  the 
frequency  of  the  pump  minus  the  frequency  of  the  weak  signal,  but  the 
pump  in  Case  II  has  a greater  frequency  than  the  pump  in  Case  I. 

The  work  presented  in  this  thesis  is  a logical  extension  of  the 
work  done  by  Schaffer  on  Case  I.  The  subject  of  this  thesis  is  the 
suppression  of  sound  by  sound  for  Case  II.  General  theoretical  models  and 
solutions  which  are  applicable  to  both  Case  I and  Case  II  are  developed  in 
Chapter  II.  The  models  and  solutions  are  compared  to  each  other  and  their 
limitations  are  discussed.  Chapter  III  is  a detailed  description  of  the 
computer  programs  which  implement  weak  shock  theory  with  provisions  for 
attenuation  and  dispersion.  The  intent  is  to  give  the  reader  enough 
insight  into  the  computer  programs  so  that  he  can  use  them  with  confidence 
in  a minimum  amount  of  time.  For  readers  who  do  not  want  or  desire  this 
insight,  Chapter  III  may  be  skipped  without  interfering  with  the  continuity 
of  the  remainder  of  the  thesis.  The  experimental  apparatus  and  its  cali- 
bration and  accuracy  are  described  in  Chapter  IV.  Theoretical  and 
experimental  results  are  presented  and  compared  in  Chapter  V.  Conclusions 
and  final  comments  are  made  in  the  last  chapter.  Chapter  VI.  Appen- 
dices A,  B,  and  C contain  listings  of  the  major  computer  programs  used. 
Appendix  D contains  a detailed  derivation  of  an  alternate  representation 
for  one  of  the  pre shock  solutions  discussed  in  Chapter  II,  the  improved 


FM  solution. 


CHAPTER  II 


r 

THEORY 

The  theoretical  tools  used  to  describe  the  interaction  between 
a finite  amplitude  tone  and  a second,  weaker  tone  are  presented  here. 
Although  the  subject  of  this  thesis  is  the  suppression  of  sound  by  sound 
for  Case  II,  much  of  the  theory  is  applicable  to  both  cases.  Whenever 
possible  a general  notation  is  employed  so  that  the  theoretical  results 
are  useful  for  both  Case  I and  Case  II.  Three  possible  approaches  were 
mentioned  in  Chapter  I to  describe  the  interaction:  Fenlon's  solution, 

an  FM  solution,  and  weak  shock  theory.  Theoretical  descriptions  of  the 
behavior  of  finite  amplitude  sound  are  frequently  limited  to  shock-free 
waveforms,  and  such  is  the  case  for  both  Fenlon's  solution  and  the  FM 
solution.  The  solution  based  on  weak  shock  theory,  however,  can  accommo- 
date waveforms  with  discontinuities.  All  three  solutions  may  be  derived 
from  an  implicit  solution  called  the  Earnshaw  solution.  In  weak  shock 
theory,  the  Earnshaw  solution  must  be  supplemented  by  a shock  propagation 
equation  derived  from  the  Rankine-Hugoniot  shock  relations.  Thus  the 
Earnshaw  solution  is  used  to  describe  the  continuous  sections  of  waveform 
between  disconu 'nuities,  and  the  shock  propagation  equation  is  used  to 
make  the  connections  across  the  discontinuities. 

In  this  chapter  the  general  form  of  the  Earnshaw  solution  is 
reviewed,  and  then  two  exact  preshock  solutions  for  a two- frequency  source 
are  given:  the  Earnshaw  solution  for  the  time  domain  and  the  Fenlon 

solution  for  the  frequency  domain.  Various  useful  approximations  of  the 
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exact  solutions  are  then  made  and  compared  with  each  other  and  with  the 
exact  solutions.  Next,  the  behavior  of  the  interaction  after  shock 
formation  is  considered.  The  shock  relations  are  combined  with  the 
Earnshaw  solution  to  form  weak  shock  theory,  and  then  the  equations  of 
weak  shock  theory  are  cast  in  difference  equation  form  and  implemented  in 
a computer  subroutine.  Next  the  importance  and  the  effect  of  attenuation 
and  dispersion  are  considered,  followed  by  mathematical  formulations  for 
attenuation  and  dispersion.  The  last  section  of  the  chapter  deals  with 
the  major  computer  algorithm,  C500,  which  calculates  the  propagation  of  a 
wave;  the  effect  of  attenuation  and  dispersion  is  included  in  the 
algorithm . 


A.  Earnshaw  Solution 

g 

The  Earnshaw  solution0  is  an  exact  solution  of  the  equations  of 
motion  for  progressive  waves.  The  solution  is  first  reviewed  for  gaseous 
fluids.  A convenient  approximation  of  the  solution  is  then  obtained  which 
is  applicable  to  liquids  as  well  as  gases.  The  equations  of  motion  for  a 

Ip 

lossless  perfect  gas  are  (see,  for  example,  Banta  ) 


where  p,  u,  and  p are  the  instantaneous  values  of  density,  particle 
velocity,  and  pressure,  respectively,  pQ  and  pQ  are  the  ambient  values  of 
density  and  pressure,  x is  propagation  distance,  t is  time,  and  y is  the 
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ratio  of  specific  heats  of  the  gas.  If  the  wave  motion  is  progressive, 
say,  restricted  to  outgoing  waves,  the  following  first  integral  of 
Eqs.  2.1  through  2.3  may  be  obtained  (for  example,  see  Blackstock15) , 


n + (c  +gu)u  =0  , 

t v O X 


(2.4) 


where  Mr+l)/2  for  gases.  A solution  of  this  equation  for  the  boundary 


condition 


u(0,t)  = f(t) 


u(x,t)  = f(cp)  , 


(2.5a) 


where 


cp  = t - 


co  + &u 


(2.5b) 


Equations  2.5a  and  2.5b  constitute  the  Earnshaw  solution.  The  parameter  cp 
represents  the  time  a point  on  a waveform,  i.e.,  a wavelet,  left  the 


origin. 


An  important  physical  interpretation  may  be  obtained  directly 
from  Eqs.  2.4  and  2.5b.  Inspection  of  these  equations  shows  that  the 


propagation  speed  is 


g = c + pu 
dt  o 


(2.6) 


In  general,  different  points,  or  wavelets,  of  a waveform  propagate  with 
different  speeds,  and  the  result  is  distortion  of  the  wave.  The  variation 


The  variation 
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in  propagation  speed  is  due  to  two  things.  First,  the  sound  speed  varies 
because  of  the  nonlinearity  of  the  pressure- density  relation;  the  sound 
speed  c is  (see,  for  example,  Blackstock  ) 


(2.7) 


where  cq  is  the  small  signal  value  of  the  sound  speed.  Second,  because 
the  propagation  of  an  acoustical  wave  is  longitudinal,  the  wave  is  con- 
vected  by  its  own  passage  through  the  medium.  The  effect  of  convection 
on  propagation  speed  may  be  expressed  as 


dx 

dt 


c + u 


(2.8) 


The  combination  of  Eqs.  2.7  and  2.8  results  in  Eq.  2.6. 

An  approximate  form  of  the  Earnshaw  solution  is  often  used  for 
convenience.  The  approximate  form  is  obtained  by  performing  a binomial 
expansion  of  Eq.  2.5b  and  keeping  terms  of  order  e only,  where  e,  the 
Mach  number,  is  defined  as  the  ratio  of  the  maximum  particle  velocity  to 
the  small  signal  speed.  The  result  of  the  binomial  expansion  is 


u = f(q>)  (2.9a) 

<P  = t'  + £2L.  f(cp)  , (2.9b) 

c 

o 

where  t'  is  the  retarded  time,  t-(x/cQ).  By  common  usage  Eqs.  2.9  are 
referred  to  as  the  Earnshaw  solution. 

The  Earnshaw  solution  is  equally  applicable  to  liquids.  A 


general  expression  for  the  equation  of  state  for  lossless  perfect  fluids 


t 
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K 

is  (for  example,  see  Beyer 


. P - po  B /p  - Po\  . C lp  ' V 

p - p = A + — r I I + TT  1 I + 

0 2 ! \ / 3-  \ P„ 


where  the  values  of  the  coefficients.  A,  B,  C,  etc.,  are  in  general 
empirically  determined.  To  the  same  order  of  approximation  as  that  used 
in  the  derivation  of  Eq.  2.9b,  the  Earnshaw  solution  is  made  applicable 
to  liquids  by  taking 

in  Eq.  2.9b.  For  air  the  ratio  B/2A  is  equal  to  (7-l)/2. 

It  is  important  to  remember  the  restrictions  on  the  validity  of 
Eqs.  2.9-  Because  of  the  approximations  made,  Eq.  2.9b  is  valid  only  for 
values  of  e«l.  Note  that  e=0.1  corresponds  to  a SPL  of  174  dB  re 
0.00002  N/m  in  air  and  a SPL  of  284  dB  re  1 pPa  in  water.  Another 
important  point  is  that  by  itself  the  Earnshaw  solution,  in  either  of  the 
forms  presented  in  Eqs.  2.5  or  2.9>  is  valid  only  for  continuous  wave- 
forms. The  solution  may  be  used  to  predict  the  formation  of  shocks. 

Beyond  the  shock  formation  distance,  however,  the  waveforms  predicted  with 
the  Earnshaw  solution  are  multivalued. 


^ 2 
Throughout  this  thesis  the  reference  pressure  for  SPL  is  = 0.00002  N/m 

for  air. 
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B.  Exact  Pre shock  Solutions  for  Two-Frequency  Source 

The  Earnshaw  solution  for  the  two- frequency  boundary  condition 


u(0,t)  = u sin  a)  t + u sin  co  t 
op  p ow  w 


u(x.t')  = u sin  cn  cp  + u sin  cjo  © 
op  p ow  wT 


(2.11a) 


cp  = t ' + / u sin  Q)  cp  + u sin  co  cp  ) . (2 . lib) 

2 \ op  p ow  w / 


The  notation  adopted  is  as  follows.  The  subscript  w denotes  the  signal 

and  the  subscript  p the  pump.  A zero  preceding  a subscript  indicates 

source  amplitude  of  the  particular  tone.  The  first  term  in  Eq.  2.11a  is 

referred  to  as  the  pump  term  u^  and  the  second  one  as  the  weak  signal 

term  u . The  combined  particle  velocity  signal,  composed  of  the  pump  and 

the  weak  signal,  is  referred  to  as  u.  When  a particular  tone  propagates 

by  itself,  that  is,  in  the  absence  of  the  other  tone,  a prime  is  used. 

For  example,  u ’ represents  the  pump  propagated  by  itself.  An  important 

frequency  parameter  is  the  ratio  fl=cu  /cu  . For  Case  I fi  is  greater  than  1, 

w p 

while  for  Case  II  fi  is  less  than  1.  The  w-p  notation  is  adopted  so  that 

the  derivations  which  follow  hold  for  both  Case  I and  Case  II. 

As  noted  in  an  earlier  section,  the  validity  of  the  Earnshaw 

solution  is  limited  to  distances  x<x,  where  x is  the  shock  formation 

q 

distance.  Fenlon  has  shown  that  the  shock  formation  distance  for  waves 
from  a two- frequency  source  is 


pe  k + fie  k 
p p w w 
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where  k=co/c  . The  region  of  validity  becomes 


or,  in  terms  of  cr. 


Pe  k 
P P 


Pe  k 

w w 


a < 1 - cr 
p w 


For  the  work  done  in  this  thesis,  a »aw.  A reasonable  approximation  of 
the  limitation  set  by  shock  formation  is 


a < 1 
P 


A convenient  way  to  display  data  as  a function  of  distance  is  to  plot  it 
against  a^.  The  theoretical  predictions  and  experimental  results  are 
displayed  in  this  fashion.  Normally  cr=l  is  the  dividing  line  between  the 
absence  and  presence  of  shocks  in  a waveform.  For  the  particular  waves 
studied  here  cr  <1  is  a good  approximation  for  the  shock-free  regions  of 
the  waveform,  but  cr  +(^<1  is  the  exact  relationship. 

The  effect  of  the  pump  on  the  weak  signal  is  expressed  by  the 
weak  signal  term  of  Eq.  2.11a, 


(u  ) = u sin 

w ow 

EARN 


co  ft' 

vl 


5K 


sin 


CD  CP 

P 


+ u 


ow 


sin 


(2.32) 


where  the  subscript  EARN  denotes  the  Earnshaw  solution.  This  portion  of 
the  Earnshaw  solution  represents  the  weak  signal  and  its  harmonics  in  the 
presence  of  the  pump.  It  also  includes  the  sidebands  resulting  from  the 
interaction  of  the  pump  with  the  weak  signal.  But,  Eq.  2.12  does  not 
describe  the  complete  sideband  formation  of  the  interaction  of  the  pump 


(B-23) 


a = ai  , 
o ’ 


y = ct/i  , 


u = y - cos  7 


In  terms  of  these  variables,  the  impulse  response  becomes 

4n2C  exp[-2a(l-y  cos  ej]  fy+1  f r ; ^ 

h(y,ro)  = 7 J — I0^2a  sin  6^1*  (y-u)  ) 

X expT-2au  cos  0^  - £Q(y2-l)/uJ  . (B-26) 

Equation  (B-26)  completes  the  reduction  of  the  impulse  response  of  a 
parametric  array  to  a single  integral.  This  integral  is  complicated,  and 
an  exact  analytic  expression  for  it  is  not  known  at  this  time.  However, 
this  formulation  of  the  impulse  response  is  useful  for  several  reasons: 

(1)  h(y,?Q)  is  always  positive  or  zero, 

(2)  h(y,rQ)  is  a smooth  function  except  at  y=l, 

(3)  the  behavior  of  h(y,?Q)  as  y approaches  1.0  and  as  y approaches 
infinity  is  well  known,  and 

(4)  the  properties  of  h(y,?  ) imply  that  it  will  be  possible  to 
approximate  h(y,?  ) by  a polynomial  in  y. 
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The  first  term  of  Eq.  2.l4  is  an  expression  for  the  amplitude  of  the  pump 
and  its  harmonics  in  the  presence  of  the  weak  signal;  the  second  term  is 
the  remaining  portion  of  u,  namely,  the  weak  signal,  its  harmonics,  and 
the  sidebands  of  the  interaction. 

The  frequencies  of  the  spectral  components  of  u are  seen  to  be 


cd  = men  + na> 
m,n  | w p 


The  indices  for  the  frequencies  are  used  to  form  a notation  for  the 

individual  frequency  components  of  the  spectra  of  u.  The  general  notation 

for  the  amplitude  of  a spectral  component  of  u is  u^  n,  where  the  angular 

frequency  of  the  component  is  equal  to  the  value  of  cd^  For  example, 

u^  q is  the  amplitude  of  the  weak  signal  component  (cd^) . For  Case  II 

(<Dp>cDw),  u^  ^ represents  the  amplitude  of  the  component  located  at  a 

frequency  of  cd  -u>w  because  a^-cD^  is  negative.  The  components  located  at 

the  frequencies  of  the  pump  and  its  harmonics  are  represented  by  uq  n 

(cd  evaluated  at  m=0  is  equal  to  nco  ) . The  first  term  of  Eq.  2.l4  is 
m,n  P 

represented  by  u and  the  second  term  by  u -u 
* J o,n  J m,n  o,n 

The  sideband  formation  is  composed  of  contributions  from  both 

terms  of  the  Earnshaw  solution,  Eq.  2.11a;  therefore,  (u  ) ^u  -u  . 

w EARN  ’ > 

Part  of  the  energy  to  form  the  sidebands  comes  from  the  suppression  of  the 
weak  signal  represented  as 


u.  _ = ——  J.(a  ) J (Act  ) 
1,0  a 1'  w o p 


(2.15a) 


= u J (fta  ) 
aw  o'  p' 


(2.15b) 


C.  Approximate  Pre shock  Solutions  for  Two-Frequency  Source 


The  approximate  result,  Eq.  2.15b,  was  used  in  Chapter  I to  predict  the 
suppression  of  the  weak  signal  for  Case  I and  Case  II  (Eq.  l.l).  Another 
source  of  energy  for  the  sideband  formation  is  the  pump  and  its  harmonics. 
The  first  summation  term  of  Eq.  2.l4  is  equal  to  J^ncr^/ft)  times  the 
Fubini  solution  for  the  pump.  (The  Fubini  solution  (for  example,  see 
Blackstock  ) describes  the  nonlinear  propagation  of  the  signal  from  a 
single  frequency  source.)  The  difference  between  the  Fubini  solution  for 
the  pump  and  uq  n of  the  Fenlon  solution  is  a measure  of  the  energy  redis- 
tributed from  the  pump  and  its  harmonics  to  the  sidebands  through  the 
interaction  of  the  pump  and  the  weak  signal.  The  difference  becomes 
greater  for  larger  n or  for  the  higher  harmonics  of  the  pump,  and  is 
larger  in  Case  II  than  in  Case  I because  fl  in  Case  II  is  less  than  1. 

For  Case  I more  energy  for  the  sideband  formation  comes  from  the  suppres- 
sion of  the  weak  signal,  whereas  in  Case  II  more  energy  is  redistributed 
from  the  harmonics  of  the  pump. 


The  exact  solutions,  Eq.  2.11  and  Eq.  2.13,  are  available  for 
analytical  investigations,  but  approximate  solutions  are  often  much  easier 
to  implement  and  just  as  useful  for  certain  purposes  as  the  exact  solu- 
tions. The  crudest  approximate  solution  is  reviewed  first  and  more 
sophisticated  approximate  solutions  are  discussed  later.  An  FM  solution 
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may  be  obtained  by  approximating  cp  (Eq.  2.11b)  in  the  Earnshaw  solution. 


If  uq^»uow^  cp  may  be  approximated  as 


cp  = t'  + u sin  cjd  cp 

c 2 op  p 

o 


(2.16) 


In  this  case  the  pump  term  of  the  Earnshaw  solution  becomes,  in  the 
frequency  domain. 


« -u  £ 


J (no  ) sin  nto  t'  , 

ncf  n p p 

P v 


(2.17a) 


which  is  recognized  as  the  classic  Fubini  solution.  In  other  words,  in 
this  approximation  the  weak  signal  has  no  effect  on  the  pump.  But  the 
pump  definitely  affects  the  weak  signal,  as  shown  by  the  expression  for  u 


u = u sin(o)  t'  + ftcr  sin  co  cp) 
w ow  w p p 


(2.17b) 


A further  simplification  is  now  made  by  letting  cp=t'  in  the  inner  sine 
function  of  Eq.  2.17.  The  result  is 


(u  ) = u sin(co  t’  + fia  sin  to  t') 


w 


FM 


ow 


(2.18) 


where  the  subscript  FM  refers  to  the  FM  solution.  Replacing  cp  by  t*  in 
this  way  is  justified  in  that  the  approximation  is  made  in  a term  that  is 
already  just  a correction  to  the  phase  of  t ' . A frequency  domain  repre- 

l8 

sentation  of  Eq.  2.18  is  (for  example,  see  Carlson  ) 


00  p 

(u  ) = u ) ' J (ftcr  ) sin|(o)  +ncn  )t'| 
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(2.19) 
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A better  approximate  solution,  which  may  be  called  the  improved 
FM  solution,  is  obtained  by  not  approximating  Cp  by  t'  in  the  inner  sine 
function  of  Eq.  2.1Jb-  The  following  Bessel  function  frequency  domain 
representation  of  Eq.  2.17b  has  been  found  by  Blackstock  (see  Appendix  D), 

(UW}T  = uow  2 fTT~n  Jn[(fi+n)ap]  sin[(VncUp)t']  ’ (2-20) 


where  the  subscript  I FM  denotes  "improved  FM  solution." 

Physical  interpretation  of  the  FM  and  improved  FM  solutions  is 
important  in  understanding  the  usefulness  of  the  solutions.  The  approxi- 
mate cp  equation  (Eq.  2.l6)  used  to  derive  these  solutions  contains  many 
implicit  assumptions:  (l)  the  pump  is  not  influenced  by  the  presence  of 

the  weak  signal  and  the  pump  term  does  not  give  rise  to  contributions  to 
the  sideband  formation  (the  pump  term  is  simply  the  Fubini  solution) , 

(2)  the  weak  signal  does  not  distort  (neither  Eq.  2.19  not  Eq.  2.20 
contains  harmonics  of  the  weak  signal),  and  (5)  the  sideband  formation  is 
due  solely  to  the  modulation  of  the  weak  signal  by  the  pump.  These 
assumptions  are  simplifying  because  it  has  been  shown^  that  both  terms  of 
the  Earnshaw  solution  contribute  to  the  sideband  formation,  the  pump  is 
affected  by  the  presence  of  the  weak  signal,  and  the  distortion  of  the 
weak  signal  is  dependent  on  its  amplitude. 

A still  better  approximate  pre shock  solution  may  be  obtained 
from  the  Fenlon  solution.  The  m=0  terms  of  Eq.  2.13b  represent  the  ampli- 
tude of  the  pump  and  its  harmonics, 


2uQp 

na  Jo 
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n > 0 


(2.21) 


If  no Jo(ncrw/U)=l  and  Eq.  2.21  reduces  to  the  Fubini  amplitudes. 
The  m=l  terms  of  Eq.  2.13b  represent  the  weak  signal  and  the  sidebands 
resulting  from  the  mutual  interaction  of  the  pump  and  the  weak  signal, 
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(2.22) 
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and  Eq.  2.22  becomes 


u..  = u J 

l,n  ow 


n[<fl+n)''p] 


(2.23) 


which  is  referred  to  as  the  approximate  Fenlon  solution  (APP  FEN) . Unlike 
the  previous  approximate  solutions,  however,  this  expression  includes 
sideband  contributions  from  both  the  pump  signal  and  the  weak  signal.  The 
m=2,3,...,  terms  of  Eq.  2.13a  represent  the  higher  harmonics  of  the  weak 
signal  and  the  sidebands  resulting  from  their  interaction  with  the  pump 
and  its  harmonics.  For  aw(m+n/ft)«l,  the  amplitudes  of  these  "higher 
order"  sidebands  are  very  small.  The  improved  FM  solution  (Eq.  2.20) 
differs  from  the  approximate  Fenlon  solution  (Eq.  2.23)  by  the  factor 


n 

fi  + n 


(2.24) 
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The  difference  between  the  two  solutions  is  small  for  Case  I,  where  H>1, 
and  for  small  n,  but  the  two  solutions  differ  considerably  for  Case  II, 
where  fl<l,  even  for  small  n.  As  already  noted,  the  reason  there  is  a 


difference  between  the  two  solutions  is  that  in  approximating  cp  by 
Eq.  2.16  the  effect  of  the  weak  signal  on  the  pump  was  not  included  in 
the  improved  FM  solution. 
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Figure  2.1  is  a flowchart  which  shows  the  relationships  between 
the  various  preshock  solutions.  In  the  figure,  equivalent  solutions  are 
across  from  one  another,  and  the  direction  of  increasing  exactness  is  from 
the  bottom  to  the  top  of  the  figure.  It  may  be  observed  from  the  figure 
that  a time  domain  solution  equivalent  to  the  approximate  Fenlon  solution 
has  not  yet  been  found. 


D.  Effect  of  the  Weak  Signal  on  the  Pump 

The  weak  signal  has  an  effect  on  the  pump,  seen  in  the  factor 
JQ(n aw/ft)  that  appears  in  the  first  summation  term  of  the  Fenlon  solution 
(Eq.  2.l4),  even  though  the  weak  wave  is  generally  regarded  as  a "small 
signal."  If  the  factor  were  unity,  the  series  would  be  precisely  the 
Fubini  solution.  An  interesting  exercise  is  to  subtract  the  Fubini  solu- 
tion, Eq.  2.17a,  that  is,  the  pump  propagated  by  itself  u^',  from  the 
complete  Fenlon  solution: 


u - u 


s 


2uqt 

na 


| /na  \ J 

lJo\  ft  / 1l 


J (na  ) 

n p 


sin 


nco  t 1 
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+ DS 


where  DS  stands  for  the  double  sum  in  Eq.  2.l4.  It  is  clear  that  u-u  ' 

P 

includes  (l)  all  the  sideband  components,  (2)  the  weak  signal  and  its 
harmonics,  and  (3)  certain  residual  signals  at  frequencies  equal  to  the 
pump  and  its  harmonics.  The  residual  signals  are  a measure  of  the  effect 
of  the  weak  wave  on  the  pump.  The  time  domain  representation  of  u-u^'  may 


i 

■ 

A 


L. 


EARNSHAW  SOLUTION  FENLON  SOLUTION 


FIGURE  2.1 

FLOWCHART  FOR  EXACT  AND  APPROXIMATE  PRESHOCK  SOLUTIONS 
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be  found  from  the  Earnshaw  solution.  The  effect  of  the  weak  signal  on  the 
pump  is  seen  in  the  time  domain  by  the  sidebands  due  to  the  pump  term  of 
the  Earnshaw  solution.  The  time  domain  representation  of  the  Fubini 
solution, 


up'(x,t')  = uQp  sin  copcp  , (2.25a) 

where 

cp  = t’  + uQp  sin  oyp  , (2.25b) 

co 

is  subtracted  from  the  complete  Earnshaw  solution  (Eq.  2.1l).  This  sub- 
tractive procedure  is  only  valid  for  the  complete  Earnshaw  solution.  In 
the  approximate  Earnshaw  solutions  the  approximation  used  for  cp  guarantees 
that  up=up';  therefore,  the  subtractive  procedure  yields  the  weak  signal 
term. 

Figure  2.2  illustrates  the  predicted  effect  of  the  pump  on  the 
weak  signal  in  Case  I using  the  following  solutions:  the  complete 

Earnshaw  solution  (Eq.  2.12),  the  FM  solution  (Eq.  2.18),  the  improved  FM 
solution  (Eq.  2.17b),  and  the  specially  constructed  signal  u-up'  (Eq.  2.11a 
minus  Eq.  2.25a),  which  shows  the  effect  of  the  weak  signal  on  the  pump. 

The  four  solutions  are,  respectively,  (a),  (b),  (c),  and  (d)  in  the  figure. 
The  source  amplitudes  and  frequencies  are  exaggerated  to  emphasize  the 
effect  of  one  signal  on  the  other. 

A comparison  of  the  spectra  in  Fig.  2.2  reveals  that  the 
predicted  amplitude  of  the  weak  tone  is  the  same  in  every  plot.  For  the 
FM  solution  the  sideband  formation  is  symmetric  about  the  weak  tone,  but 
for  the  improved  FM  solution  the  formation  is  asymmetric.  The  higher 
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frequency  sidebands  are  enhanced.  In  the  complete  Earnshaw  solution  the 
weak  signal  undergoes  finite  amplitude  distortion  at  a rate  depending  on 
its  amplitude,  but  in  the  EM  solutions  the  weak  signal  does  not  distort  no 
matter  what  its  amplitude.  The  greater  high  frequency  sideband  formation 
in  the  Earnshaw  solutions  is  partially  a result  of  the  harmonics  of  the 
weak  signal  interacting  with  the  pump  and  its  harmonics.  The  difference 
between  plot  (a)  and  plot  (b)  is  more  difficult  to  explain.  Both  the  pump 
and  weak  signal  terms  in  the  Earnshaw  solution  contribute  to  the  sideband 
formation  caused  by  the  mutual  interaction  between  the  tones.  In  plot  (a), 
only  the  weak  signal  term  of  the  Earnshaw  solution  is  shown;  therefore  only 
the  sidebands  due  to  the  effect  of  the  pump  on  the  weak  signal  are  seen. 

The  Earnshaw  solution  in  plot  (d)  shows  the  effect  of  the  pump  on  the 
weak  signal  and  the  effect  of  the  weak  signal  on  the  pump.  The  observed 
sideband  formation  is  due  to  both  terms  in  the  solution.  The  fundamental 
and  harmonics  of  the  pump  (the  residual  signals  described  earlier)  seen  in 
the  plot  indicate  that  the  presence  of  the  weak  signal  has  an  effect  on 
the  finite  amplitude  distortion  of  the  pump,  as  is  very  evident  in  the 
waveform  of  plot  (d).  The  weak  signal  seems  to  be  undergoing  amplitude  as 
well  as  frequency  modulation.  The  amplitude  modulation  is  not  caused  by 
the  relatively  large  amplitude  of  the  weak  signal  in  this  example;  compu- 
tations with  weak  signals  of  much  smaller  amplitude  (100  dB  down  from  pump 
amplitude)  showed  similar  amplitude  modulation. 

Figure  2.3  illustrates  the  predicted  effect  of  the  pump  on  the 
weak  signal  in  Case  II.  Plot  (a)  illustrates  the  Earnshaw  solution 
(Eq.  2.12),  plot  (b)  the  FM  solution  (Eq.  2.18),  plot  (c)  the  improved  FM 


t'«t-x/c0 

Ml 

00  - 

1 

00  H 

1 

| 

0 

4 ' 6 

12  16 

FREQUENCY  - kHi 

o)  EARNSHAW  SOLUTION 

, Up  REMOVED 

t « t-x/e„  j 


0 4 S 12  16 

FREQUENCY  - kHi 

W F.M.  SOLUTION,  up  REMOVED 


/ ♦'•♦-x/c0  J< 

V) 

AO  - 

00  - 

0 

4 1 

12  16 

FREQUENCY  . kH. 

(e)  IMPROVED  F.M.  SOLUTION,  up  REMOVED 

» ■ t-x/c.  jJ 


0 4 I 12  16 

FREQUENCY  - kHi 

(<•)  EARNSHAW  SOLUTION,  u'p  REMOVED 

FIGURE  2.3 

EFFECT  OF  PUMP  ON  WEAK  SIGNAL  CASE  H 
SOURCE  SPL„  * 135  4B,  f w • 1 kHi 
SOURCE  SPLp  * 150  4B,  fp  * 10  kHi 
X * 0.7  m 


35 


solution  (Eq.  2.17b),  and  plot  (d)  the  signal  u-u  The  source 

P 

amplitudes  and  frequencies  are  exaggerated  to  emphasize  the  interaction 
between  the  tones. 

Inspection  of  the  spectra  of  Fig.  2.3  reveals  that  the  predicted 
amplitude  for  the  weak  tone  is  the  same  for  all  the  solutions.  The  FM 
solution  sidebands  are  symmetric  about  the  pump  and  harmonic  frequencies; 
the  improved  FM  solution  sidebands  are  asymmetric  (by  a fraction  of  a 
decibel) . In  the  improved  FM  solution  the  sidebands  around  the  second 
harmonic  of  the  pump  are  20  dB  larger  than  those  for  the  simple  FM 
solution.  The  finite  amplitude  distortion  of  both  tones  is  incorporated 
in  the  complete  Eamshaw  solution  but  not  in  the  approximate  solutions. 

The  distortion  of  the  weak  signal  is  observed  in  the  Earnshav  solution 
plots  [plots  (a)  and  ( d) ] - The  interaction  of  the  pump  with  the  harmonics 
of  the  weak  signal  is  a cause  of  the  greater  sideband  formation  in  the 
Earnshaw  solutions.  The  sideband  formation  in  plot  (a)  is  due  solely  to 
the  weak  signal  term  of  the  two- frequency  Earnshaw  solution.  The  sidebands 
in  plot  (d)  are  due  to  both  terms  of  the  Earnshaw  solution.  Evidence  of 
the  effect  of  the  weak  signal  on  the  pump  is  the  energy  located  at  fre- 
quencies equal  to  that  of  the  pump  and  its  second  harmonic.  Notice  that 
the  amplitude  of  the  second  harmonic  component  is  larger  than  that  of  the 
pump.  The  sidebands  in  plot  (d)  are  abnormally  large  compared  to  those  in 
the  other  solutions.  Only  the  sidebands  from  the  weak  signal  term  of  the 
Earnshaw  solution  are  contained  in  the  other  solutions.  The  effect  of 
both  tones  on  each  other  is  best  observed  in  the  dramatic  waveforms  of 


plot  (d).  The  period  of  the  cusps  is  approximately  the  period  of  the 
pump  (or  the  sidebands  observed  in  the  corresponding  spectrum).  The  cusp 
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formation  is  not  due  to  the  relatively  large  amplitude  of  the  weak  signal 

in  Fig.  2.3;  cusps  were  predicted  in  calculations  where  the  amplitude  of 

the  weak  signal  was  100  dB  less  than  the  pump  amplitude. 

Figures  2.2  and  2.3  demonstrate  that  the  primary  source  of 

energy  for  the  sidebands  is  different  in  the  two  cases.  The  suppression 

of  the  weak  signal  for  both  cases  is  described  by  Eq.  2.15a.  The 

difference  between  the  Fubini  solution  and  the  u terms  of  the  Fenlon 

o,n 

solution  (Eq.  2.14-)  is  a measure  of  the  amount  of  energy  of  the  pump  and 
its  harmonics  which  is  redistributed  to  the  sidebands.  The  normalized 
difference  between  the  two  is 
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°,n  FEN  °,n  FUB 
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FUB 


(2.26) 


where  the  subscript  FUB  refers  to  the  Fubini  solution,  n=l  corresponds  to 
the  pump,  and  n=2  to  the  pump's  second  harmonic,  etc.  The  smaller  the 
difference  between  the  solutions,  the  less  pump  energy  is  redistributed  to 
the  sidebands.  For  Case  I,  the  argument  of  Jq  in  Eq.  2.15a  may  be  greater 
than  one  (a^<l  but  ft>l)  and  may  equal  the  first  zero  of  Jq.  Therefore, 
much  energy  may  be  redistributed  from  the  weak  signal  to  the  sidebands. 

On  the  other  hand,  for  Case  I little  energy  is  redistributed  from  the  pump 
and  its  harmonics  so  long  as  n is  not  large;  Eq.  2.26  is  approximately 
zero  because  Jo(na^/(l)±l  for  ctw«1  and  0>1.  For  Case  II  the  roles  of  the 
primary  and  secondary  sources  of  energy  for  the  sideband  formation  are 
switched.  The  argument  of  Jq  in  Eq.  2.15  has  a limiting  value  of  1 
(fKl  and  c?p<l)  for  Case  II,  so  that  little  suppression  of  the  weak  signal 
is  achieved.  On  the  other  hand,  the  amount  of  energy  redistributed  from 
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the  pump  and  its  harmonics  is  greater.  Even  though  a «1,  the  fact  that 

w 

ft<l  means  that  the  argument  of  JQ  in  Eq.  2.26  may  deviate  substantially 
from  zero  and  thus  the  value  of  JQ  may  deviate  correspondingly  from  unity. 

The  energy  for  the  sideband  formation  in  both  Case  I and  Case  II 
comes  primarily  from  the  higher  frequency  tone  and  its  harmonics.  In 
Case  I the  higher  frequency  tone  is  the  weak  signal  and  so  a great  deal 
of  its  energy  is  required  to  form  the  sidebands.  In  Case  II  the  pump  is 
the  higher  frequency  tone.  The  sideband  energy  taken  from  it  represents 
such  a small  percentage  that  the  relative  decrease  in  the  amplitude  of 
the  pump  and  its  harmonics  is  not  very  great. 

Figure  2.4  illustrates  the  difference  between  the  approximate 
Fenlon  solution  and  the  Bessel  function  representation  of  the  improved  FM 
solution  due  to  Blackstock  (see  Appendix  D) . The  approximate  Fenlon 
solution  spectra  represent  u^  n and  the  FM  solution  represents  u^.  The 
ratio  of  the  improved  FM  solution  to  the  approximate  Fenlon  solution 
(Eq.  2.24)  represents  the  difference  between  the  two  solutions.  Plots  (a) 
and  (b)  show  that  for  Case  I there  is  little  difference  between  the  two 
solutions  (the  ratio  is  approximately  one  for  small  n and  ft>l),  but  for 
Case  II  the  difference  is  larger  (with  fl<l,  even  for  small  n,  the  ratio  is 
not  equal  to  one).  For  Case  II  the  predicted  sidebands  for  the  approximate 
Fenlon  solution  are  much  larger  than  those  for  the  improved  FM  solution. 
Similarly  large  sidebands  were  also  observed  in  Fig.  2.3(d).  The  two 
solutions  with  large  sideband  formations  both  incorporate  to  some  degree 
the  mutual  sideband  formation  of  the  pump  and  the  weak  signal. 


SPL  - dB  5-  SPL 
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FREQUENCY  - kHz  FREQUENCY  - kHz 


APPROXIMATE  FENLON  SOLUTION,  u,,  „ (b)  IMPROVED  FM  SOLUTION,  uw 

CASE  I 

SOURCE  SPLp  = 150  dB,  fp  = 1 kHz 
SOURCE  SPLW  = 135  dB,  fw  = 10  kHz 
x - 1.7  m 


FREQUENCY  - kHz  FREQUENCY  - kHz 


(c)  APPROXIMATE  FENLON  SOLUTION,  u,  „ 

'#  n 


(d)  IMPROVED  FM  SOLUTION,  uw 


CASE  n 

SOURCE  SPLW  s 135  dB,  fw  = 1 kHz 
SOURCE  SPLp  = 150  dB,  fp  = 10  kHz 
x = 0.7  m 


L.  ^ 


FIGURE  2.4 

COMPARISON  OF  FENLON  AND  IMPROVED  FM  SOLUTIONS  FOR 
THE  EFFECT  OF  PUMP  ON  WEAK  SIGNAL,  CASE  I AND  CASE  II 
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E.  Shock  Relations 

When  shocks  form  in  a waveform,  conservation  equations  must  be 
applied  across  them  to  connect  the  continuous  portions  of  the  waveform  on 
either  side.  The  conservation  equations  commonly  used  are  called  the 
Rankine-Hugoniot  shock  relations.  The  relations  may  be  combined  to  yield 
an  expression  for  the  propagation  speed  of  a shock.  The  expression  which 
is  an  approximation  of  the  same  order  as  the  Earnshaw  solution  is  (for 

g 

example,  see  Pestorius  ) 


shock 


(Ug+Ui) 


(2.27) 


where  ^-1  is  the  propagation  speed  of  the  shock,  the  subscript  2 

/shock 

refers  to  behind  the  shock,  and  the  subscript  1 refers  to  in  front  of  the 


shock.  The  propagation  velocity  of  a wavelet  in  a continuous  waveform  is 

given  by  Eq.  2.6.  A comparison  of  the  two  equations  (Eqs.  2.6  and  2.27) 

shows  that  a shock  travels  with  the  mean  speed  of  wavelets  just  behind  and 

just  ahead  of  it.  Therefore  wavelets  behind  the  shock  tend  to  catch  up 

with  the  shock  while  the  shock  tends  to  overtake  wavelets  ahead  of  it. 

The  Earnshaw  solution  is  valid  only  for  progressive  waves.  Is 

this  condition  satisfied  when  a wavelet  overtakes  a shock,  that  is,  is 

there  a reflection  of  the  wavelet  which  would  negate  the  progressive  wave 

19 

assumption?  Lighthill  calculated  that  for  weak  shocks  the  departure 
from  progressive  wave  flow  is  of  the  order  where  Z is  the  shock  strength 

O 

defined  as  (p^-p^/p^.  Thus,  to  order  Z the  formation  of  shocks  does  not 
invalidate  the  progressive  wave  assumption.  When  e=0.1,  the  maximum  value 


of  Z is  0.0032.  Therefore,  for  the  purpose  of  this  thesis,  reflections 
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from  shocks  may  be  ignored.  The  Earnshaw  solution  may  continue  to  be 
used  for  the  continuous  sections  of  a waveform. 

If,  when  a wavelet  overtakes  a shock,  or  vice  versa,  there  is  no 
reflection,  what  is  the  consequence  of  the  merger?  The  answer  to  this 
question  is  an  important  result  of  nonlinear  acoustics.  A shock  is  com- 
posed of  at  least  two  wavelets.  One  of  the  wavelets  is  associated  with 
the  peak  (maximum  pressure)  of  the  shock  and  another  wavelet  is  associated 
with  the  trough  (minimum  pressure)  of  the  shock.  These  two  wavelets, 
which  may  be  used  to  describe  a shock,  may  be  referred  to  as  the  endpoint 
wavelets  of  the  shock.  In  a wavelet-shock  merger  either  the  shock  or  the 
wavelet  overtakes  the  other  because  of  differing  propagation  speeds.  A 
wavelet-shock  merger  is  really  the  merging  of  one  of  the  shock's  endpoint 
wavelets  with  the  third  wavelet.  When  a wavelet  overtakes  a shock,  the 
wavelet  governs  the  amplitude  of  the  shock,  but  when  a shock  overtakes  a 
wavelet  ahead  of  it  the  overtaken  wavelet  determines  the  amplitude  of  the 
shock.  A similar  process  occurs  between  two  shocks  with  differing  propa- 
gation speeds  which  overtake  one  another.  The  result  of  the  merging  of 
two  shocks  is  a single  shock. 

F.  Weak  Shock  Theory 

Weak  shock  theory  is  based  upon  the  Earnshaw  solution  (Eq.  2.9) 

and  the  equation  for  the  propagation  speed  of  a shock  (Eq.  2.27) • The 

quantity  cp  in  the  Earnshaw  solution  has  units  of  time,  and  the  coefficient 
2 

3f(cp)/co  has  units  of  time  per  distance.  The  shock  propagation  relation- 
ship connects  the  Earnshaw  solution  on  one  side  of  a shock  to  the  Earnshaw 
solution  on  the  other  side  of  the  shock.  It  would  be  very  useful  if  the 
shock  propagation  equation  were  in  units  of  time  per  distance.  Inverting 
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u 


Eq.  2.27  will  produce  the  desired  dimensions.  Further  usefulness  is 
achieved  by  making  a variable  transformation  from  t,  real  time,  to  t ' , the 
retarded  time.  The  variable  transformation  is  expressed  by 


dt'  = dt  _ _1 
dx  dx  c 


(2.28) 


Substituting  the  inverse  of  Eq.  2.29  in  this  equation  yields 


shock 


_1_ 

c 

o 


(2.29) 


Performing  a binomial  expansion  to  the  same  order  as  before  and  collecting 
terms  yields 


dt 

dx 


) 


shock 


~^~~2  (Ul+U2} 
2c 

o 


(2.30) 


Equations  2.30  and  2-9  form  weak  shock  theory. 

In  preparation  for  the  implementation  of  weak  shock  theory  in  a 
computer  algorithm,  the  following  ideas  need  to  be  established.  A wave- 
form may  be  described  as  a collection  of  wavelets  associated  with  a 


corresponding  collection  of  retarded  times,  or  as 


u1(t1  ) , i — 1,  2,  3,  ..,  n (2.31) 

where  n is  the  number  of  points  in  the  waveform  representation.  To 
simplify  this  discussion  and  other  descriptions,  let  time  t refer  to 
retarded  time  t ' . The  use  of  weak  shock  theory  to  propagate  the  wave 
represented  by  Eq.  2.31  results  in  a new  waveform  represented  by  the  same 


, 


1*0 

collection  of  wavelets,  but  a different  collection  of  times.  The  new 
collection  of  times  is  calculated  from  Eqs.  2.9b  and  2. 30.  The  changed 
time  relationships  between  the  wavelets  represents  the  distortion  of  the 
waveform.  A shock  is  formed  when  two  or  more  wavelets  are  associated  with 
the  same  time  or  when  a multivalued  region  is  predicted  in  the  waveform. 

G.  Weak  Shock  Computer  Algorithm 

Equations  2.30  and  2.9b  form  the  mathematical  basis  for  weak 
shock  theory.  The  application  of  the  equations  was  briefly  discussed  in 
the  last  section,  but  how  is  weak  shock  theory  practically  implemented? 
Analytical  implementation  is  possible  but  usually  very  difficult  to  do 
except  for  the  simplest  of  boundary  conditions.  In  principle,  implemen- 
tation of  weak  shock  theory  by  means  of  a computer  algorithm  is  simple. 

Q 

Pestorius  has  developed  a computer  algorithm  referred  to  as  WAVEPROP 
which  implements  weak  shock  theory.  WAVEPROP  was  initially  designed  to  be 
used  on  random  waveforms,  but  works  equally  as  well  on  deterministic 
wave  forms . 

In  WAVEPROP  an  input  waveform  is  propagated  a desired  distance 
in  incremental  steps.  The  input  waveform  is  represented  as  a collection 
of  wavelets  and  their  corresponding  times.  The  values  of  the  wavelets  are 
stored  in  an  array,  called  the  U array,  and  the  values  of  t'  in  a second 
array,  called  the  T array.  In  the  propagation  of  the  initial  waveform  the 
first  incremental  distance  produces  a new  T array  associated  with  the  U 
array.  The  spacing  between  the  elements  of  the  T array  is  initially 
uniform  but,  because  of  waveform  distortion,  the  spacing  between  the  T 
elements  changes  with  every  incremental  propagation  step.  The  input  for 
the  next  propagation  step  is  the  output  of  the  last  step.  A shock  is 

_ 1 


formed  when  two  or  more  elements  of  the  U array  are  associated  with  the 
same  time.  If  more  than  two  wavelets  art  associated  with  the  same  value 
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of  time,  only  the  largest  and  smallest  amplitude  wavelets  (the  endpoint 
wavelets)  are  needed  to  represent  the  shock.  The  remaining  wavelets 
associated  with  the  shock  are  not  necessary  and  are  dropped  from  the 
arrays.  Therefore,  as  a waveform  undergoes  propagation,  the  number  of 
elements  in  the  U and  T arrays  decreases. 

The  form  of  the  weak  shock  theory  equations  implemented  by  the  1 

computer  programs  is  not  Eqs.  2.30  and  2.9b.  Because  of  the  step-by-step 
nature  of  the  propagation  process,  Eq.  2. 9b  may  be  viewed  in  a new  light. 

The  elements  of  the  T array  before  a propagation  step  are  referred  to  as 
t'  ^ and  the  elements  of  the  array  after  a propagation  step  are  referred 
to  as  t1  With  these  definitions  and  a little  rearranging,  Eq.  2.9b 

becomes 


t' 

new 


t'  old 


ef(t'oid)^ 


(2.32) 


where  f(t'  is  the  value  of  the  wavelet  associated  with  t'  _ , and  Ax  is 

the  size  of  the  incremental  propagation  step.  Weak  shock  theory  in  the 
form  of  Eqs.  2.30  and  2.32  is  in  a convenient  form  to  be  implemented  in 
the  computer  algorithm  discussed.  All  that  is  needed  is  to  express 

3 

Eqs.  2.30  and  2.32  in  difference  equation  form  (due  to  Pestorius  ).  The 
difference  equation  representation  of  Eq.  2.32  is 


- 
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and  the  difference  equation  representation  of  Eq.  2.30  is 

UgJkAx,  tg ' (kAx) | + u^jkAx,  tg'(kAx) 

(2. 3*0 

where 

Ax  = incremental  distance, 

t'[(k+l)Ax]  = the  value  of  t'  after  k+1  incremental  distance  steps, 

t = the  value  of  t1  associated  with  a shock,  and 
s 

t ’[(k+l)Ax]  = the  value  of  t ' after  k+1  incremental  distance  steps, 
s s r 

Equation  2-33  may  he  rearranged  and  substituted  into  Eq.  2.34  to  obtain  a 
simplified  shock  propagation  equation 

t '[(k+l)Ax]  + t _'[(k+l)Ax] 

tg'[(k+l)Ax]  g- • (2.35) 

Equations  2.33  and  2.35  are  the  form  of  weak  shock  theory  used  by  WAVEPROP 
to  calculate  the  new  times  for  a collection  of  wavelets. 

H.  Attenuation  and  Dispersion 

Attenuation  and  dispersion  need  to  be  incorporated  in  the 
computer  algorithms  because  the  predictions  generated  are  to  be  compared 
to  experimental  results.  The  experiments  are  performed  in  a plane  wave 
tube,  and  the  presence  of  the  tube  wall  introduces  appreciable  attenuation 
and  dispersion.  In  a normal  sized  duct,  or  a pipe,  the  effect  of 
viscosity  and  heat  conduction  near  the  bounding  surfaces  is  the  dominant 
cause  of  attenuation  over  a wide  frequency  range.  In  a pipe  there  exists 
a viscous  boundary  layer  across  which  the  particle  velocity  associated 
with  an  acoustical  signal  undergoes  a rapid  transition  from  its  mainstream 


PAxc 

t *[(k+l)Ax]  = t '(kAx) — 

S S d. 


h 
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value  to  zero  at  the  boundary.  A thermal  boundary  also  exists  across 
which  the  propagation  of  an  acoustical  signal  changes  from  adiabatic  in 
the  mainstream  to  isothermal  at  the  boundary.  These  boundary  layers  are 
usually  very  thin.  The  mainstream  is  the  part  of  the  interior  of  the  pipe 
where  the  boundary  layer  effects  are  not  felt  by  the  propagating  acoustic 
signal.  Figure  2.5  is  a plot  of  boundary  layer  thickness  versus  frequency 

for  air  for  both  the  thermal  boundary  and  the  viscous  boundary.  The 
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viscous  boundary  layer  (Y^)  is  given  by 


Y 

v 


(2.56) 


where  v is  the  kinematic  viscosity  of  the  medium.  The  thermal  boundary 
layer  thickness  (Y  ) is  equal  to  Y divided  by  the  square  root  of  the 

w V 

Prandtl  number.  In  1868  Kirchhoff  published  his  derivation  of  an  expres- 
sion for  the  mainstream  attenuation  coefficient,  a , and  one  for  the 

’ nr 
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boundary  layer  attenuation, 


“b  = 


1 + 


y - 1 

y/pr 


rc 


vco 

2 


(2.57) 


o 


where  Pr  is  the  Prandtl  number  defined  as  c u/k,  \ is  the  dilatational 
viscosity  coefficient,  p is  the  coefficient  of  sheer  viscosity,  r is  the 
radius  of  the  pipe,  k is  the  coefficient  of  thermal  conductivity,  and  c^ 

is  the  specific  heat  at  constant  pressure.  If  (^>>0:^  and  r»Yv,  a pipe  is 

22 

referred  to  by  Weston  as  a wide  pipe.  For  the  plane  wave  tube  used 


BOUNDARY  LAYER  THICKNESS 
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(5  cm  diam)  and  for  a tone  of  5 kHz,  the  ratio  of  o^/a^  is  equal  to  0.005. 
The  next  most  important  source  of  attenuation  is  relaxation  effects  in  the 
medium.  Attenuation  due  to  relaxation  is  equal  to  the  attenuation  due  to 
boundary  layer  effects  (afe)  at  frequencies  of  the  order  of  100  kHz.  ^ 
Therefore,  for  the  plane  wave  tube  used  and  at  the  frequencies  of  the 
primaries  (less  than  5 kHz),  the  major  cause  of  attenuation  of  the  acoustic 
signals  is  viscous  boundary  layer  effects. 

The  viscous  and  thermal  boundary  layers  cause  dispersion  as  well 
as  attenuation.  In  air  for  tubes,  dispersion  is  expressed  by 


c(co)  = 


(2.59) 


For  a simple  frequency  source,  the  various  harmonics  get  out  of  step  with 
each  other  because  of  their  varying  phase  velocities.  Only  for  large  fre- 
quencies is  the  phase  velocity  equal  to  cq.  The  net  result  is  that  the 

peak  of  a shock  becomes  rounded  while  the  trough  stays  sharp.  An  illus- 

24 

tration  due  to  Webster  of  the  effect  of  dispersion  on  a sawtooth  is 
shown  in  Fig.  2.6. 


I.  Major  Computer  Algorithm 


Program  C500,  the  computer  algorithm,  which  implements  weak 
shock  theory  with  attenuation  and  dispersion,  is  fairly  involved.  A 
detailed  explanation  of  program  C500  is  presented  in  Chapter  III,  and  a 
documented  listing  of  the  program  and  the  subroutines  used  are  presented 
in  Appendix  A.  A brief  outline  is  given  here  so  the  reader  may  gain  a 


basic  understanding  of  the  program. 


Figure  2.7  is  a simplified  flowchart  for  program  C 500.  The 


program  begins  with  the  initialization  of  program  parameters  including 
reading  from  data  cards,  the  source  SPLs,  and  frequencies  of  the  two  source 
tones  and  the  desired  propagation  distance.  The  initial  waveform  is 
represented  by  a collection  of  wavelets,  whose  values  are  stored  in  the 
U array,  and  a corresponding  collection  of  times,  which  are  stored  in  the 
T array.  The  calculation  of  the  initial  waveform  takes  place  next,  and 
then  the  waveform  is  stored  on  disk.  The  initial  waveform  is  propagated 
the  desired  distance  by  a series  of  incremental  propagation  steps. 

DO  loop  80l6  is  the  major  DO  loop  of  the  program;  each  completion  of  the 
loop  signifies  the  completion  of  one  incremental  propagation  step.  DO 
loop  8016  begins  by  recalling  the  waveform  (the  U and  T arrays)  from  disk 
storage  and  writing  the  arrays  on  magnetic  tape  to  be  displayed  later  by 
program  TIMEPLOT.  The  propagation  of  the  waveform  is  achieved  by  calling 
subroutine  WAVEPROP,  the  weak  shock  computer  algorithm.  The  output  wave- 
form of  WAVEPROP  is  then  stored  on  disk.  A call  is  made  to  subroutine 
TIMEFRE,  which  performs  a fast  Fourier  transform  (FFT)  on  the  waveform  so 
that  the  Fourier  coefficients  of  the  waveform  may  be  calculated.  The 
largest  Fourier  coefficients  are  then  selected  and  their  amplitudes  and 
frequencies  corrected  for  inaccuracies  inherent  in  the  FFT.  Next,  the 
corrected  Fourier  coefficients  are  written  on  a second  magnetic  tape  to  be 
processed  later  by  program  GRAPH.  In  the  next  section  of  the  program  a 
decision  is  made  about  applying  attenuation  and  dispersion  corrections. 

The  corrections  are  made  only  every  NATDIS  number  of  trips  through 
DO  loop  8016.  If  attenuation  and  dispersion  corrections  are  not  made. 


control  of  the  program  returns  to  the  beginning  of  DO  loop  8016.  The  next 


PROGRAM  C500 


FIGURE  2.7 

SIMPLIFIED  FLOWCHART  FOR  PROGRAM  C500, 
THE  WEAK  SHOCK  COMPUTER  ALGORITHM 
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and  last  section  of  the  program  applies  any  corrections  that  are  to  be 
made.  Here  the  waveform  is  recalled  from  disk  storage  and  an  FFT 
obtained  (by  a call  to  subroutine  TIMEFKE)  because  the  corrections  are 
made  in  the  frequency  domain.  An  inverse  FFT  takes  the  signal  back  into 
the  time  domain.  The  corrected  waveform  is  stored  on  the  disk  storage  and 
control  of  the  program  returns  to  the  beginning  of  DO  loop  80l6.  The 
whole  process  in  DO  loop  8016  is  repeated  until  the  desired  propagation 
distance  is  reached.  An  important  aspect  of  program  0500  is  the  use  of 
disk  storage.  When  attenuation  and  dispersion  corrections  are  not  made, 
the  disk  storage  allows  subroutine  WAVEPROP  to  have  the  output  of  the  last 
propagation  as  its  input  for  the  next  propagation  step.  Example  program 
outputs  are  presented  in  Chapter  III. 

Program  C500  is  used  to  calculate  the  theoretical  predictions 
for  comparison  to  the  experimental  results  in  Chapter  V.  In  Chapter  III 
the  output  of  program  C500  is  compared  to  one  of  the  exact  pre shock 
solutions  for  the  same  source  conditions.  The  comparison  (Fig.  3-6)  shows 
that  the  results  of  the  two  solutions  are  equivalent;  therefore,  nothing 
is  lost  by  using  program  C500  for  x<x^.  As  a matter  of  fact,  the  use  of 
program  C500  is  advantageous  because  it  provides  a simple  method  of 
incorporating  attenuation  and  dispersion. 


CHAPTER  III 


MAJOR  COMPUTER  ALGORITHMS 

The  computer  algorithms  used  to  implement  weak  shock  theory  with 
provisions  for  attenuation  and  dispersion  are  explained  in  detail  here. 
This  chapter  is  intended  to  help  a U6er  of  the  programs  become  familiar 
with  the  scope  and  limitations  of  the  programs.  A reader  not  desiring 
this  insight  may  skip  the  present  chapter  without  disrupting  the  con- 
tinuity of  the  thesis.  Program  C500,  which  implements  the  weak  shock 
method  and  also  includes  provisions  for  attenuation  and  dispersion,  is 
discussed  first.  After  the  detailed  flowchart  for  the  program  is 
explained,  a discussion  of  the  major  subroutines  and  problems  associated 
with  the  program  is  given.  Also  programs  GRAPH  and  TIMEPLOT,  which 
display  the  output  of  program  C500,  are  discussed  and  their  flowcharts 
explained;  samples  of  the  output  of  programs  GRAPH  and  TIMEPLOT  are 
included.  The  calibration  of  the  three  programs  as  a unit  is  discussed 
next.  In  the  last  section  of  the  chapter  the  flexibility  of  these 
programs  is  considered  for  applications  other  than  the  very  precise  ones 
dealt  with  in  this  thesis. 

A.  Major  Computer  Algorithms 

Flowchart  for  Program  C^OO 

It  may  be  useful  to  recall  the  simplified  flowchart  for 
program  C500  found  in  Chapter  II  (Fig.  2.4)  before  discussing  the  detailed 
flowchart  for  program  C500,  Fig.  3.1.  That  simplified  flowchart  and 


accompanying  discussion  serve  as  a good  introduction  to  program  C500  and 
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FIGURE  3.1 

DETAILED  FLOWCHART  FOR  PROGRAM  C500 


U(I+1) 


FIGURE  3.1  (Con t.) 

DETAILED  FLOWCHART  FOR  PROGRAM  C500 
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this  chapter.  Another  aid  in  understanding  the  flowchart  for  program  C500 
is  Appendix  A,  which  is  a documented  listing  of  program  C500  and  its 
subroutines  and  definitions  for  the  various  program  parameters. 

The  program  parameters  that  determine  the  output  of  program  C500 
are  read  from  data  cards  at  the  beginning  of  the  program.  The  definitions 
of  the  parameters  are  found  in  Appendix  A.  An  important  relationship 
between  DISTANCE  (the  desired  propagation  distance)  and  DISINC  (the  incre- 
mental propagation  6tep  size)  is  that  the  ratio  DISTANCE/DISINC  must  be 
equal  to  an  integer.  Next,  the  source  SPL  of  each  tone  (called  A or  B in 
the  program)  is  used  to  obtain  a source  particle  velocity  amplitude.  The 
vast  majority  of  program  parameters  are  then  initialized.  Some  explana- 
tion, besides  the  definitions  in  Appendix  A of  a few  parameters,  is 
needed  and  will  help  to  explain  some  basic  idea  behind  the  program.  The 
number  of  points  in  the  initial  waveform  is  the  even  number  NPOINTS.  The 
initial  waveform  includes  an  even  number  NPD  of  pump  wavelengths  to  help 
ensure  zero  or  approximately  zero  endpoints.  NDISTOR  is  the  number  of 
incremental  propagation  steps  done  in  succession  without  stopping  to 
Fourier  transform  the  output  waveform.  NO  is  the  number  of  (NDISTORS)x 
(DISINC)  propagation  steps  necessary  to  reach  the  desired  propagation 
distance.  The  initial  waveform  is  represented  by  a particle  velocity 
array  U and  a time  array  T.  The  initial  waveform  is  calculated  by 
DO  loop  100.  The  phase  of  the  initial  waveform  is  l80°  to  ensure  that  the 
shocks  which  form  will  form  inward  away  from  the  endpoints  of  the  arrays. 
Zero  buffers  are  added  to  the  beginning  of  the  waveform  (DO  loop  101)  and 
to  the  end  of  the  waveform  (DO  loop  102)  to  allow  for  possible  spreading 
out  of  the  original  time  base  caused  by  distortion  when  the  wave 


propagates.  The  initial  waveform,  which  represents  the  boundary- 
condition,  is  stored  on  disk  so  that  a waveform  may  be  recalled  and 
operated  on  (e.g.,  transformed  by  FFT)  without  destroying  the  original 
waveform.  Some  of  the  input  (deterministic)  parameters  and  other  program 
parameters  are  written  on  the  spectrum  output  magnetic  tape  (magnetic 
tape  No.  l)  to  identify  the  data  on  the  tape. 

The  next  section  of  the  program  consists  of  the  remainder  of  the 
program  and  contains  the  major  DO  loop  of  the  program,  DO  loop  8oi6.  An 
incremental  propagation  of  the  waveform  stored  on  disk  is  achieved  by  the 
completion  of  one  trip  through  DO  loop  80l6;  every  completion  of  a loop 
represents  the  completion  of  a propagation  step  equal  to  (NDISTORS)x 
(DISINC) . NO  number  of  loops  of  DO  loop  80l6  have  to  be  completed  in 
order  to  propagate  an  initial  waveform  the  desired  distance.  The  DO  loop 
has  three  major  sections:  the  first  calculates  the  Fourier  transform  of 

the  output  waveform  of  a propagation  step;  in  the  second,  the  major  (local 
maxima)  Fourier  coefficients  are  selected  from  the  output  waveform  spec- 
trum and  corrected  values  and  frequencies  are  calculated  for  the  selected 
coefficients;  and  the  last  section  contains  the  implementation  of  the 
attenuation  and  dispersion  corrections.  The  DO  loop  begins  by  incrementing 
DISTANCE  by  (DISINC)* (NDISTOR)  for  the  calculation  of  (referred  to  as 
SIG  in  the  program)  for  the  next  output  waveform.  The  waveform  on  disk 
is  recalled  and  written  on  magnetic  tape  No.  2 to  be  plotted  by  program 
TIMEPLOT  later,  if  desired.  A call  to  subroutine  WAVEPROP  achieves  the 
propagation  of  the  waveform  through  the  weak  shock  method.  The  output 
waveform  returned  by  WAVEPROP  is  then  stored  on  disk.  A fast  Fourier 
transform  (FFT)  is  performed  on  the  output  waveform  by  a call  to 
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subroutine  TIMEFKE . Subroutine  TIMEFRE  removes  the  zero  buffers  from  the 

waveform  and  calls  subroutine  RESAMPLE.  Subroutine  RESAMPLE  resamples  the 

waveform  to  produce  a desired  number  of  equally  spaced  points  in  the  U and 

T arrays.  Still  under  the  control  of  subroutine  TIMEFRE,  the  resampled 

waveform  undergoes  an  FFT  and  the  real  and  imaginary  parts  of  the  Fourier 

coefficients  are  returned  to  the  main  program.  A data  window  is  applied 

to  the  frequency  components  of  the  waveform  by  a call  to  subroutine 

HANGEN2.  A data  window  is  used  to  ensure  that  the  endpoints  of  the  wave- 

* 

form  are  zero  to  reduce  side  lobes,  aliasing,  and  leakage  inherent  in  a 
FFT.  The  real  and  imaginary  parts  of  the  frequency  components  are  scaled 
to  calibrated  values  by  DO  loop  20.  Calibrated  values  are  necessary  so 
that  the  output  of  program  C500  is  compatible  with  experimental  results. 

DO  loop  315  calculates  the  magnitude  of  the  real  and  imaginary  points  of 
the  frequency  components  and  DO  loop  888  calculates  their  frequencies. 

The  sections  of  the  program  between  DO  loop  888  and  the  end  of  DO  loop  315 
compute  the  first  and  last  Fourier  coefficient  of  the  output  waveform. 

In  the  next  section  of  the  program,  and  of  DO  loop  8016,  the 
major  or  local  maximum  frequency  components  of  the  spectrum  of  the  wave- 
form are  selected.  The  chosen  frequency  components  are  corrected  for  the 
data  window  and  for  the  inherent  error  of  the  FFT.  The  corrected  ampli- 
tude and  frequency  of  a Fourier  coefficient  are  calculated  by  an  interpo- 

25 

lation  scheme  by  Burgess.  The  selection  of  the  maximum  frequency 
components  is  accomplished  by  DO  loop  5>  and  is  initiated  by  checking  the 
frequency  of  a coefficient.  If  the  frequency  of  the  coefficient  is 

greater  than  FMAX  (the  largest  frequency  of  interest),  DO  loop  5 is 

_ 

To  be  defined  and  discussed  in  the  FFT  subroutine  section. 
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terminated  and  the  remaining  section  of  the  program  is  executed.  When  a 
maximum  frequency  component  is  located,  the  interpolation  is  between  the 
maximum  component  and  the  larger  of  the  two  components  to  the  left  and 
right  of  the  maximum.  The  corrected  amplitude  of  a frequency  component 
is  stored  in  an  even  numbered  element  of  the  FT  array  and  the  frequency 
of  the  component  is  stored  in  the  next  larger  odd  element  of  the  FT  array. 
The  array  containing  the  corrected  Fourier  coefficients  and  frequencies  is 
written  on  magnetic  tape  No.  1,  with  the  value  of  associated  with  the 
spectrum.  Magnetic  tape  No.  1,  or  the  spectrum  output  tape,  is  processed 
by  program  GRAPH. 

The  corrections  for  attenuation  and  dispersion  are  applied  in 

t the  last  section  of  the  program  and  the  third  section  of  DO  loop  80l6. 

The  corrections  are  not  applied  after  every  propagation  step  because  to 
do  so  would  introduce  errors  due  to  resampling  (resampling  problems  are 
discussed  in  the  subroutine  RESAMPLE  section).  When  the  corrections  are 
not  applied,  the  control  of  the  program  execution  is  returned  to  the 
beginning  of  DO  loop  8016.  An  important  point  about  program  C500  is  that 
when  attenuation  and  dispersion  corrections  are  not  applied,  the  input 
waveform  for  the  next  propagation  step  is  the  output  of  the  last  propaga- 
tion step  which  has  been  stored  on  disk.  DO  loop  815  implements  the 
corrections  for  attenuation  and  dispersion  only  when  a given  number, 
called  NATDIS,  of  completions  of  DO  loop  8016  have  been  executed.  Before 
the  corrections  are  made,  the  waveform  is  recalled  from  disk  storage  and 
subroutine  TIMEFRE  is  called  to  transform  the  signal  to  the  frequency 
domain.  The  corrections  are  made  and  then  an  inverse  FFT  is  taken  to  get 
back  into  the  time  domain.  DO  loop  8l6  scales  the  corrected  waveform  to 
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calibrated  values.  The  waveform  amplitude  has  to  be  calibrated  so  that 
the  finite  amplitude  distortion  predicted  by  weak  shock  theory  for  the 
wave  is  correct.  DO  loop  8l8  shifts  the  elements  of  the  U and  T arrays 
forward  to  make  room  for  the  right-hand  zero  buffer  which  was  removed  by 
subroutine  TIMEFRE.  A zero  buffer  is  added  to  the  right  end  of  the 
corrected  waveform  by  DO  loop  821;  a zero  buffer  at  the  left  end  is  added 
by  DO  loop  817.  Next,  the  corrected  waveform  with  zero  buffers  is  stored 
on  disk,  and  the  signal  is  returned  to  the  beginning  of  DO  loop  8016  until 
the  desired  propagation  distance  is  reached. 

a.  Subroutine  WAVEPROP 

Subroutine  WAVEPROP  is  the  weak  shock  method  subroutine.  An 
input  waveform,  described  in  terms  of  a U array  and  a T array,  is  propa- 
gated an  incremental  distance  according  to  weak  shock  theory.  WAVEPROP 
renumbers  the  elements  of  the  arrays  as  shock-wavelet  merging  causes  a 
decrease  in  the  total  number  of  elements  in  the  arrays.  Equally  spaced 
data  points  are  not  required.  The  original  version  of  WAVEPROP  performed 
all  the  necessary  DISINC  steps  to  propagate  the  initial  waveform  the 
desired  distance  with  one  call  to  WAVEPROP.  The  present  form  of  WAVEPROP 
performs  NDISTORS  number  of  DISINC  steps  with  every  call  to  WAVEPROP.  For 

g 

a complete  explanation  of  WAVEPROP  see  Pestorius. 

b.  Subroutine  TIMEFRE 

Subroutine  TIMEFRE  was  specifically  written  for  C500.  TIMEFRE 
performs  the  actual  transformation  from  the  time  domain  to  the  frequency 
domain.  The  zero  buffers  are  checked  to  see  if  they  are  large  enough. 

If  they  are,  the  zero  buffers  are  removed  from  the  U array  prior  to  a call 
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to  the  FFT  and  the  elements  of  the  U and  T arrays  are  renumbered,  but  the 
time  relationships  between  the  elements  of  the  arrays  are  maintained.  If 
the  zero  buffers  are  not  large  enough  (they  no  longer  exist),  the  program 
execution  is  stopped  and  one  of  the  following  messages  is  printed: 

RIGHT  ZERO  BUFFER  NOT  LARGE  ENOUGH 
if  the  right-hand  zero  buffer  is  not  large  enough,  or 
LEFT  ZERO  BUFFER  NOT  LARGE  ENOUGH 

if  the  left-hand  zero  buffer  is  too  small.  It  is  then  necessary  to 
increase  N00  or  NOOL  depending  on  which  of  the  zero  buffers  is  in  question. 
N00  and  NOOL  are  the  number  of  zeroes,  respectively,  in  the  right  and  left 
zero  buffers.  The  number  of  remaining  zeroes  in  the  two  buffers  is 
printed  with  every  call  to  TIMEFRE.  NET  is  the  number  of  zeroes  in  the 
left  buffer,  and  LET  in  the  right  buffer.  TIMEFRE' s calling  sequence  is 
CALL  TIMEFRE  (U,  T,  NPOINTS,  N,  NEXP2,  IX,  SCAL,  NH, 

ITOT,  TMAX,  LU,  II,  SIG,  MAGIF,  DT), 

where 

N is  the  desired  number  of  elements  in  U or  T after  a call 
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LU  is  a look-up  table  for  the  FET, 

II  is  the  number  of  the  last  incremental  step  taken, 

MAGIF  is  a dummy  variable, 

DT  is  the  spacing  in  the  time  data  after  RESAMPLE, 
and  the  remaining  variables  have  been  previously  defined.  A call  to 
TIMEFKE  causes  a call  to  RESAMPLE  and  a call  to  IFFTDF.  Subroutine 
TIMEFRE  returns  the  real  and  imaginary  parts  of  the  Fourier  coefficients 
to  the  main  program  in  the  U array. 

c.  Subroutine  RESAMPLE 

Subroutine  RESAMPLE  used  in  program  C500  is  very  nearly 

g 

identical  to  the  version  used  by  Pestorius.  The  only  difference  is  that 
the  internal  storage  for  the  computational  steps  has  been  changed  to  allow 
the  storage  area  to  be  used  by  other  parts  of  program  C500.  Besides 
respacing  the  data  points,  calls  to  RESAMPLE  have  an  unfortunate  side 
effect.  In  the  process  of  obtaining  more  (or  less)  array  elements  which 
are  equally  spaced,  the  waveform  defined  by  the  arrays  becomes 
"undi storted. " The  distortion  of  a waveform  is  represented  by  the  changed 
spacing  of  the  elements  of  the  T array:  the  elements  near  a shock  become 
bunched  up  and  elements  near  a rarefaction  are  spread  apart.  A call  to 
RESAMPLE,  however,  produces  even  spacing  in  the  T array;  even  with  an 
interpolation  scheme,  the  resulting  waveform  is  undistorted  (the  shocks 
are  not  as  sharp).  In  program  C500  the  undistortion  of  a waveform  is  a 
problem  only  when  attenuation  and  dispersion  are  applied.  When  attenua- 
tion and  dispersion  are  not  applied,  the  input  to  WAVEFROP  is  the  previous 
output,  which  has  been  stored  on  disk.  To  minimize  the  RESAMPLE  problem, 
not  every  Incremental  propagation  step  is  followed  by  attenuation  and 
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dispersion  corrections.  Pestorius  has  a detailed  discussion  on  selecting 
a value  of  NATDIS,  the  number  of  DISINC  between  corrections  for  attenua- 
tion and  dispersion. 


d.  FET  Subroutine 

The  Fourier  transform  routine,  the  forward  transform  IFFTDF  or 
the  backward  transform  IFFTDB,  is  an  integer  routine.  The  reason  an 
integer  FET  is  used  on  a floating  point  array  is  that  it  saves  computa- 
tional time.  Before  a call  to  IFFTDF  is  made  in  subroutine  TIMEFRE,  the 
U array  is  scaled  up  by  SCAL  so  as  to  maintain  the  desired  number  of 
significant  digits  before  the  truncation  involved  in  the  transition  from 
floating  point  to  integer  number  occurs.  The  E'ET  routines  are  designed 
for  harmonic  analysis  of  the  frequency  composition  of  waveforms.  The 
frequencies  of  interest  in  a two-frequency  interaction  are  in  general  not 
harmonically  related.  An  FET  routine  may  be  used  to  find  nonharmonically 
related  frequency  information  by  making  the  fundamental  sought  by  the  FET 
small.  The  hope  is  that  the  harmonics  of  the  fundamental  will  fall  near 
the  true  component  frequencies.  Using  a small  fundamental  is  a "shotgun" 
approach  that  may  lead  to  erroneous  results.  If  a true  component  fre- 
quency does  not  fall  directly  in  the  middle  of  the  main  lobe  of  a harmonic 
of  the  fundamental,  the  amplitude  of  the  Fourier  coefficient  returned  by 
the  FET  will  be  incorrect.  Instead  of  the  correct  value  for  the  frequency 
component,  two  false  components  are  returned  by  the  FET.  A sin(x)/x 
interpolation  is  performed  between  the  two  false  values  to  get  the  correct 

amplitude  and  frequency  of  the  component.  The  interpolation  scheme  used 

23 

also  corrected  for  the  Hanning  data  window.  If  the  length  of  a given 
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set  of  time  data  is  T,  the  fundamental  frequency  Af  associated  with  an  FFT 
of  the  data  is 

Af  = l/T 

A long  length  of  data  is  required  for  a small  Af. 

Another  problem  encountered  with  an  FFT  is  aliasing.  An  FFT  has 
a maximum  frequency  for  which  a Fourier  coefficient  can  be  correctly 
calculated.  The  frequency  of  the  last  Fourier  coefficient  that  may  be 
calculated  with  an  FFT  is  called  the  Nyquist  or  cutoff  frequency. 

| Nyquist  Frequency  = l/(2At)  , 

where  At  is  that  sampling  rate  given  by  At=T/n,  where  n is  the  number  of 
points  in  the  time  data.  If  the  time  data  has  real  frequency  components 
which  are  greater  than  the  Nyquist  frequency,  aliasing  occurs.  A fre- 
quency greater  than  the  Nyquist  frequency  is  folded  over  to  a frequency 
less  than  the  Nyquist  frequency.  If  the  Nyquist  frequency  is  10  kHz  and 
there  is  a real  frequency  component  at  15  kHz,  the  Fourier  transform  will 
fold  the  15  kHz  into  a 5 kHz  component  in  the  spectrum  of  the  waveform. 

The  folding  over  of  frequency  components  greater  than  the  Nyquist  frequency 
by  a Fourier  transform  may  be  represented  by 

F = F - (M)*(F  ) 

apparent  true  ' v Nyquist'' 

where  F is  the  folded  over  frequency,  F^.rue  is  the  real  frequency, 

^Nyquist  *8  Nyquist  frequency  for  the  Fourier  transform  for  a certain 
input,  and  M is  a positive  integer  such  that  ^’apparen-t  Is  between  0 Hz  and 
the  Nyquist  frequency.  Aliasing  is  a potential  problem  in  the  present 
work  because  of  the  large  amplitudes  of  the  pump  harmonics;  the  more 
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remote  harmonics  may  exceed  the  Nyquist  frequency.  However,  the  problem 
is  really  important  only  when  the  folded  over  frequency  falls  upon  a fre- 
quency of  interest.  Aliasing  may  be  detected  by  changing  the  Nyquist 
frequency  by  a known  amount;  this  change  will  cause  any  aliasing  to  be 
shifted  a corresponding  amount. 

Leakage  is  also  a problem  with  an  FFT.  Leakage  is  caused  when 
the  values  of  particle  velocity  associated  with  endpoints  of  the  time  data 
are  not  the  same,  causing  false  frequency  components  to  appear  throughout 
the  spectrum  of  the  waveform.  The  result  of  an  FFT  of  a finite  length 
waveform  is  not  a true  Fourier  transform  of  the  waveform  but  a transform 
of  a waveform  in  which  the  original  waveform  is  periodically  extended  to 
infinity  in  the  time  domain.  Discontinuities  at  the  periodically  extended 
waveform  cause  leakage,  and  leakage  causes  side  lobes  to  form  around  the 
main  Fourier  coefficients.  One  way  to  solve  the  leakage  problem  is  to 
ensure  that  the  signal  is  zero  at  each  endpoint.  The  initial  time  data  is 
constructed  so  that  there  is  an  integral  number  of  periods  of  the  pump. 

The  phase  of  the  input  data  is  selected  so  that  the  shocks  which  are 
formed  by  the  pump  form  inward,  away  from  the  endpoints  of  the  signal. 

When  a dispersion  correction  is  applied  to  an  intermediate  waveform,  the 
phase  of  the  frequency  components  of  the  waveform  is  changed.  The 
endpoints  of  the  waveform  are  then  not  necessarily  zero.  Zero  buffers  are 
then  required  to  allow  the  nonzero  endpoints  of  the  signal  to  distort  out 
of  the  time  base  (beyond  the  time  endpoints).  A Hanning  data  window  is 
applied  to  the  data  to  ensure  zero  endpoints  before  an  FFT  is  taken.  But 
data  windows  distort  the  data  they  are  applied  to;  therefore,  the  data  have 
to  be  corrected  after  a data  window  has  been  used.  With  all  the  problems 
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inherent  with  a FFT,  it  would  seem  to  be  a very  poor  analytical  tool. 

Burgess  has  put  it  nicely:  "It  is  the  interpretation  of  a spectral 

2*5 

representation  that  is  important,  not  its  appearance."  J 
e.  Subroutine  HANGEN2 

Subroutine  HANGEN2  applies  a Hanning  data  window  to  a waveform, 
so  as  to  shape  the  data  to  give  a desired  characteristic.  A common  data 
window  is  the  boxcar  window  which  is  applied  by  multiplying  a finite 
length  of  time  domain  data  by  a constant  function  of  magnitude  one.  The 
effect  of  applying  a boxcar  window  on  a length  of  data  results  in  simply 
the  original  length  of  data.  A Hanning  data  window  ensures  that  the  end- 
points of  a finite  length  of  data  are  zero  and  that  the  transition  to  zero 
at  the  endpoints  of  the  data  is  gradual.  A Hanning  window  is  achieved  by 
multiplying  the  time  domain  data  with  the  first  half  of  a cycle  of  a sine 
function  squared.  The  same  result  may  be  obtained  through  convolution  in 
the  frequency  domain  of  the  transformed  time  domain  data  and  the  sine 
squared  data  window;  HANGEN2  uses  the  convolution  approach.  Subroutine 
HANGEN2  and  IFFTDF  are  library  routines  of  the  Computer  Science  Division 
of  Applied  Research  Laboratories,  The  University  of  Texas  at  Austin, 

Box  8029,  Austin,  Texas  78712. 

2.  Flowchart  for  Program  GRAPH 

The  spectra  of  the  waveforms  produced  by  the  incremental 
propagation  of  an  initial  waveform  in  program  C500  are  processed  and 
displayed  in  an  appropriate  format  by  program  GRAPH.  Figure  3-2  is  a 
detailed  flowchart  for  program  GRAPH.  Important  program  parameters  from 
program  C500  which  identify  the  data  on  the  magnetic  tape  are  first  read 


FIGURE  3.2 
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from  the  program  C500  output  spectrum  tape;  program  GRAPH  variables  are 
initialized  in  the  next  section  of  the  program  (Appendix  B is  a documented 
listing  of  program  GRAPH  which  includes  the  definitions  of  the  program 
variables).  A few  of  the  parameters  warrant  additional  explanation.  A 
numerical  filter  is  employed  to  select  specific  frequency  components  from 
the  spectra  of  program  C500;  XYZ  is  the  bandwidth  of  the  numerical  filter. 
The  amplitude  for  a particular  frequency  component  is  found  as  a function 
of  . Each  selected  frequency  component  is  associated  with  a spectral 
number:  the  lowest  frequency  component  has  a spectral  number  of  1,  the 

next  lowest  2,  etc.  The  output  of  GRAPH  consists  of  plots  of  the  ampli- 
tude of  the  frequency  components  corresponding  to  the  spectral  numbers  as 
a function  of  cr^.  DO  loop  6005  calculates  the  center  frequencies  for  the 
numerical  filter  and  stores  them  in  the  W array;  the  center  frequencies 
are  really  the  frequencies  of  interest  in  the  interaction  of  the  weak 
signal  and  the  higher  frequency  pump  (o^o^a^tno^;  n an  integer).  The 
numerical  filtering  is  accomplished  by  DO  loop  200,  DO  loop  100  being  the 
actual  filter.  Each  trip  through  loop  200  is  started  by  reading  a spec- 
trum from  the  magnetic  tape.  Each  spectrum  is  searched  for  frequency 
components  which  fall  within  XYZ  Hz  from  the  center  frequencies.  If  more 
than  one  frequency  component  is  found  within  XYZ  Hz  from  a center  fre- 


quency, the  component  which  has  the  greatest  amplitude  is  chosen.  The 
selected  frequency  components  are  stored  in  the  AMP  and  SIGMA  arrays 
according  to  their  spectral  numbers.  For  example,  the  amplitude  of  the 
third  Fourier  coefficient  found  within  the  bandwidth  of  the  second 
smallest  center  frequency  is  stored  in  the  AMP(2,3)  element  of  the  AMP 
array,  and  the  corresponding  value  of  for  this  coefficient  is  stored 
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in  the  SIGMA(2,3)  element  of  the  SIGMA  array.  Each  time  a component  is 
associated  with  a spectral  number,  an  index  for  that  spectral  number  is 
incremented  by  one,  the  final  value  of  an  index  equaling  the  number  of 
selected  components  for  that  particular  spectral  number.  The  indexes  are 
called  Ml,  JG,  M3,  and  M4  for  the  first  four  spectral  numbers.  The 
remaining  portion  of  program  GRAPH,  DO  loop  801,  does  the  plotting.  The 
program  parameters  LPST  and  LPEND  determine  the  number  of  plots  to  be 
created  by  program  GRAPH.  The  first  plot  corresponds  to  a spectral  number 
equal  to  LPST,  and  the  last  plot  has  a spectral  number  equal  to  LPEND. 

The  actual  plotting  of  arrays  AMP  and  SIGMA  is  straightforward 
and  needs  no  explanation.  However,  before  the  elements  of  array  AMP  for  a 
particular  center  frequency  are  plotted,  a polynomial  least  squares  best 
fit  of  the  data  is  performed.  The  parameter  KK  in  program  GRAPH  is  equal 
to  the  order  of  the  polynomial  fit;  KK  had  a value  of  10  here.  A poly- 
nomial fit  of  the  data  is  necessary  for  two  reasons.  Attenuation  and 
dispersion  corrections  are  not  applied  at  every  incremental  propagation 
step,  but  when  corrections  are  made,  a stair-step  decrease  is  caused  in 
the  amplitudes  of  all  the  frequency  components.  The  polynomial  fitting 


The  second  reason  a polynomial  fit  is  necessary  is  that  it  connects 
elements  of  an  array  corresponding  to  a frequency  component  in  a natural, 
reasonable  manner.  Every  spectrum  searched  will  not  produce  frequency 
components  for  every  spectral  number;  therefore,  ]a~ger  jumps  might  exist 
in  between  consecutive  elements  in  SIGMA  than  caused  by  the  propagation 
increment.  Nothing  is  known  between  two  consecutive  points  in  the  AMP 


array  except  the  points  are  expected  to  be  joined  by  a smooth  curve. 

Figure  3*3  shows  an  example  output  from  program  GRAPH. 

3.  Flowchart  for  Program  TIMEPLOT 

The  waveforms  produced  by  the  incremental  propagation  of  an 
initial  waveform  by  WAVEPROP  in  program  C500  are  displayed  in  an  appro- 
priate format  by  program  TIMEPLOT.  Program  TIMEPLOT  simply  removes  the 
left  zero  buffers  from  selected  waveforms  and  plots  the  first  portion  of 
the  particle  velocity  waveform  as  a function  of  time.  Figure  J.k  is  a 
flowchart  for  program  TIMEPLOT;  Appendix  C gives  a documented  listing  of 
the  program,  including  definitions  for  the  program  parameters.  The 
program  parameters  are  initialized  in  the  beginning  of  the  program.  Two 
of  these  parameters,  IHOP  and  TMAX,  require  a little  explanation.  IHOP-1 
is  equal  to  the  number  of  waveforms  skipped  between  the  waveforms  which 
are  plotted,  and  TMAX  is  the  length  of  the  portion  of  the  plotted  waveform. 
Loop  100  is  the  major  loop  of  program  TIMEPLOT.  The  loop  is  begun  by 
reading  a waveform  from  a magnetic  tape.  The  index  of  loop  100  is  com- 
pared to  IHOP,  and  if  the  two  are  equal  the  waveform  is  plotted  and  IHOP 
is  increased  by  an  amount  equal  to  its  original  value;  if  the  two  are  not 
equal  the  index  of  loop  100  is  incremented  and  the  next  waveform  is  read 
from  the  magnetic  tape.  After  a particular  waveform  has  been  selected  to 
be  plotted,  the  zeroes  in  the  left  zero  buffer  are  counted  by  loop  308- 
The  final  value  of  the  loop  index  JJ,  less  one,  is  equal  to  the  number  of 
zeroes  in  the  left  zero  buffer.  DO  loop  310  removes  the  left  zero  buffer; 
TSLIP  is  the  amount  the  time  axis  is  slipped  backward  when  the  left  zero 
buffer  is  removed  to  retain  the  correct  time  relationship  between  the 
elements  of  the  remaining  waveform.  DO  loop  555  determines  the  number  of 


PROGRAM  TIMEPLOT 


MHOP  = IHOP  = 3 
TMAX  - (.002) (.66) 
TSCAL  = TMAX  15. 

TTIC  « .00002 
AMP  « 4. 

AMPSCL  = 2. (AMP/2.5) 
HOR  = 0 
1FLTT  = 1 


II  = l,NO 


JJ  = 1 ,N 


I = 1,NN 


MA  = I+JJ-1 
U(I)  = U (MA) 

T(I)  - T(MA)  - TSLIP 


AXES  ARE  DRAWN  AND  TIME  DOMAIN  DATA  IS  PLOTTED 


FIGURE  3.4 

FLOWCHART  FOR  PROGRAM  TIMEPLOT 
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elements  of  the  waveforms  necessary  to  plot  a length  of  the  waveform  equal 
to  TMAX.  Then  the  selected  portion  of  the  first  part  of  the  waveform  is 
plotted,  and  control  of  the  program  returns  to  the  beginning  of  loop  100. 

Figure  3.5  is  an  example  output  from  program  TIMEPLOT. 

B.  Calibration  of  Computer  Programs 

The  absolute  calibration  of  the  programs  is  necessary  because 
the  results  of  the  programs  are  to  be  compared  to  experimental  results. 

The  calibration  of  C500  concerned  the  forward  and  backward  FFTs  taken 
during  the  execution  of  C500.  An  initial  waveform  composed  of  a tone  of 
known  SPL  was  propagated  a distance  equal  to  zero.  The  zero  propagation 
distance  meant  that  the  waveform  would  experience  no  distortion  or 
attenuation  but  would  be  exposed  to  the  vast  majority  of  mathematical 
procedures  a signal  experiences  when  it  is  run  through  the  program.  The 
amplitude  of  the  Fourier  coefficient  for  the  input  tone  from  the  spectrum 
of  the  output  of  C500  was  compared  to  the  input  amplitude.  The  frequency 
of  the  tone  was  changed  so  that  the  entire  range  of  interest  in  the  fre- 
quency domain  could  be  checked.  For  a wide  range  of  input  amplitudes, 
the  input  and  output  amplitudes  for  the  tone  differed  by  no  more  than 
0.01  dB.  The  next  stage  of  the  calibration  procedure  was  to  input  a 
waveform  consisting  of  two  tones  whose  amplitude  and  frequency  were  simi- 
lar to  those  used  in  a routine  experiment  for  Case  II.  The  difference  in 
input  and  output  amplitudes  for  zero  propagation  distance  was  always  found 
to  be  within  0.02  dB.  The  amplitudes  of  the  Fourier  coefficients  corre- 
sponding to  the  sideband  frequencies  located  at  the  sum  and  difference 
frequencies  were  at  least  70  dB  down  from  the  pump  when  the  pump  was 
149  dB  and  the  propagation  distance  was  zero.  The  final  stage  of  the 
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calibration  was  to  compare  the  results  of  program  C500  to  the  results  of 
an  analytical  solution  for  the  same  input  conditions.  Figure  3-6  is  a 
comparison  of  the  predicted  suppression  for  Case  II  (before  shock  forma- 
tion) for  program  C500  and  the  FM  solution,  both  without  losses.  The 
agreement  between  the  two  is  excellent — within  0.02  dB. 


C.  Flexibility  of  Computer  Programs 

Computer  program  C500  may  be  easily  used  for  purposes  other  than 
the  specific  ones  in  this  thesis.  The  program  was  designed  to  calculate 
the  propagation  of  a Case  II  waveform.  The  propagation  of  a Case  I wave- 
form may  be  calculated  by  the  program  if  the  pump  source  SPL  and  frequency 
are  input  into  the  programs  as  B and  12,  respectively,  and  the  weak  tone 
source  SPL  and  frequency  are  input  as  A and  FI,  respectively.  The  propa- 
gation of  a pure  tone  may  be  calculated  if  the  frequency  and  source  SPL 
of  the  tone  are  input  as  F2  and  B,  respectively,  and  A is  assigned  a value 
of  zero.  If  corrections  for  attenuation  and  dispersion  are  not  desired, 
they  may  be  bypassed  by  adding  the  statement  "GO  TO  8016"  to  the  program, 
two  lines  beyond  statement  701. 

Program  C500  may  also  be  used  to  compute  the  propagation  of  more 
complicated  signals,  such  as  noise.  The  size  of  the  zero  buffers  needed 
depends  on  the  particular  signal  to  be  propagated.  For  random  noise  the 
length  of  the  two  buffers  may  have  to  be  the  same  as  the  length  of  the 
data.  An  arbitrary  waveform  could  be  introduced  in  the  program  and  stored 
in  the  U and  T arrays  in  any  number  of  ways:  data  cards,  an  analytical 

expression,  digitized  analog  data,  etc.  Program  C500  is  flexible  in  what 
constitutes  an  acceptable  input  waveform;  program  GRAPH  is  flexible  in  the 
frequency  components  it  searches  for.  The  center  frequencies  of  the 
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FIGURE  3.6 

COMPARISON  OF  PREDICTIONS  FOR  THE  SUPPRESSION  OF 
THE  WEAK  SIGNAL  FOR  CASE  H (NO  LOSSES  INCLUDED) 
SOURCE  SPLW-  115  dB,  fw  ■ 1 kHz 
SOURCE  SPLp  * 150  dB,  fp  = 4.5  kHz 
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numerical  filter  in  program  GRAPH  are  stored  in  the  W array.  The 

amplitude  of  any  desired  frequency  component  may  he  plotted  as  a function 

of  by  changing  the  values  stored  in  the  elements  of  the  W array.  For 

a single  or  double  frequency,  waveform  a is  calculated  using  the  corre- 

P 

sponding  e of  the  source  SPL  input  as  B.  For  noise  or  other  complicated 
signals,  cr^  is  not  meaningful;  the  frequency  components  of  these  signals 
may  be  plotted  as  a function  of  distance  by  redefining  SIG  in  program 
C500  as  DISTANCE.  If  individual  spectra  are  required,  the  spectra  which 
are  written  on  magnetic  tape  may  be  printed  or  plotted  easily. 
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CHAPTER  IV 

EXPERIMENTAL  APPARATUS 

The  experimental  apparatus  and  procedures  are  described  here, 
and  the  experimental  measurements  are  reported  in  the  next  chapter.  The 
experimental  apparatus  is  divided  into  three  parts:  the  transmit  system, 

the  plane  wave  pipe,  and  the  receive  system.  The  equipment  used  in  each 
system  is  listed  and,  if  necessary,  its  function  explained.  Special 
problems  experienced  with  individual  pieces  of  equipment  are  discussed. 

The  second  half  of  the  chapter  includes  the  checking  and  testing  of  each 

* 26 
part  of  the  experimental  apparatus.  Figure  4.1  (adapted  from  Webster  ) 

is  a general  diagram  for  the  experimental  setup. 

A.  The  Systems 

1.  The  Transmit  System 

The  following  list  of  equipment  was  used  in  the  transmit  system: 

a.  Bruel  and  Kjaer  (B&K)  type  1022  beat  frequency 
oscillator  (BFO), 

b.  General  Radio  (GR)  type  1310-A  oscillator, 

c.  Dukane  type  1A921  200  W power  amplifiers, 

d.  University  type  ID-75  75  W horn  driver, 

e.  Hewlett  Packard  53OOB  timing  system, 

f.  Hewlett  Packard  350C  attenuator  set, 

g.  Tektronix  51+5B  dual  channel  oscilloscope, 

h.  JBL  Model  375-H  horn  drivers  (with  aluminum  diaphragm) , 

i.  JBL  Model  375“H  horn  driver  (with  phenolic  diaphragm), 
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j.  Pearson  Electronics,  Inc.,  model  110  current  loop, 

k.  Hewlett  Packard  400  EL  PET  voltmeter,  and 

l.  Hewlett  Packard  3580A  spectrum  analyzer. 

The  transmit  system  was  further  divided  into  two  separate  subsystems: 
one  to  produce  the  pump  signal  and  another  to  produce  the  weak  signal. 
Separate  oscillators,  amplifiers,  and  drivers  were  used  to  produce  the  two 
signals  independently  to  decrease  the  intermodulation  distortion  (IM) 
between  the  signals  in  the  transmit  system.  Figure  4.2  is  a diagram  of 
the  transmit  system  showing  the  use  of  the  listed  equipment.  Some  comment 
about  the  two  types  of  JBL  horn  drivers  is  necessary.  The  JBL  375-H 
driver  with  an  aluminum  diaphragm  has  a flat  frequency  response  from 
300  Hz  to  8600  Hz,  with  a maximum  input  of  no  more  than  25  W.  The 
JBL  375-H  driver  with  a phenolic  diaphragm  has  a flat  frequency  response 
from  300  Hz  to  4.5  kHz,  with  a maximum  input  of  60  W.  Two  different 
arrangements  of  drivers  were  used  to  produce  the  two  acoustic  signals. 

In  one  case  the  weak  signal  was  produced  by  a University  ID-75  driver  and 
the  pump  signal  by  the  phenolic  diaphragm  JBL  375-H  driver.  In  the  second 
arrangement  an  aluminum  diaphragm  JBL  375-H  driver  produced  the  weak 
signal,  and  three  drivers  of  the  same  type,  operating  in  parallel,  pro- 
duced the  pump  signal.  The  first  arrangement  of  drivers  was  used  in  the 
first  experiments.  After  the  failure  of  the  phenolic  diaphragm  JBL  375-H 
driver,  the  second  arrangement  of  drivers  was  used. 

2.  Plane  Wave  Pipe 

The  plane  wave  tube  was  used  to  ensure  planar  progressive 

propagation  of  the  acoustic  signals.  This  particular  tube  has  been  used 

6 8 24  27 

in  several  previous  investigations.  ’ ’ ’ The  tube  consists  of  eight 
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3-66  m (12  ft)  lengths  of  aluminum  pipe  with  an  inside  diameter  of  5 cm. 
Both  ends  of  each  pipe  have  specially  made  connecting  flanges  affixed  with 
three  set  screws.  Each  set  of  flanges  has  three  steel  alignment  pins 
which  help  ensure  a smooth  junction  between  two  consecutive  pieces  of 
pipe.  Each  pair  of  connecting  flanges  also  contains  a microphone  port. 
Figure  4-. 3 is  a picture  of  one  of  the  flange  assemblies.  When  a micro- 
phone port  is  not  in  use,  an  aluminum  plug  is  inserted  in  the  port.  When 
a microphone  is  to  be  used  in  a certain  port,  a teflon  collar  is  inserted 
in  the  port.  A B&K  l/4  in.  microphone  may  then  be  inserted  in  the  collar 
so  that  the  microphone  diaphragm  is  flush  with  the  inside  surface  of  the 
pipe.  The  eighth  section  of  pipe  has  a 2 m tapered  fiberglass  wedge.  The 
termination  helps  satisfy  the  progressive  wave  approximation. 

The  two  different  arrangements  of  drivers  used  to  produce  the 
acoustical  signal  required  two  different  driver  mounts  for  the  plane  wave 
tube.  When  a two-driver  source  was  used,  a single  Y yoke  was  employed  to 
connect  the  drivers  to  the  end  flange.  A double  Y yoke  was  used  when  a 
four-driver  source  was  required.  The  yokes  made  a gradual  transition  from 
the  dimensions  of  the  mouth  of  each  driver  to  the  dimensions  of  the  inside 
of  the  pipe.  The  propagation  distances  of  the  observed  waveforms  are 
measured  from  the  point  in  the  yokes  at  which  the  complete  pump  signal  is 
formed.  The  single  yoke  is  shown  in  Fig.  4.4a,  and  the  double  yoke  in 
Fig.  4.4b.  In  both  pictures  the  first  microphone  location  (called  the 
monitor  microphone  location),  with  a l/4  in.  microphone  inserted,  may  be 
seen.  In  neither  picture  are  the  drivers  connected  to  the  power  ampli- 
fiers, but  the  black  conductor  in  the  lower  left  corner  is  the  trans- 
mission line  from  the  amplifiers.  Both  pictures  also  show  that  the  entire 
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FIGURE  4.3 

EXPLODED  VIEW  OF  FLANGE  ASSEMBLY  WITH 
MEASURING  MICROPHONE  AND  PLUG 

(THE  PROTECTIVE  GRID  SHOWN  ON  THE  MICROPHONE 
WAS  REMOVED  BEFORE  THE  MICROPHONE  WAS 
INSERTED  IN  THE  HOLDER  AND  FLANGE) 


(Taken  from  Schaffer*^) 
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plane  wave  tube  is  suspended  from  the  ceiling  in  a hall.  Noise  control 
treatment  of  the  plane  wave  tube  was  necessary  because  of  leakage  of  the 
acoustic  signals.  In  Fig.  4.4  the  noise  control  treatment  may  be  seen 
(fiberglass  insulation  and  duct  seal). 

The  separate  electronic  transmit  systems  did  not  completely 
eliminate  undesired  mixing  of  the  two  source  signals.  The  drivers 
connected  through  acoustical  paths  were  acting  as  receivers  of  the  acous- 
tical  signals  of  the  other  drivers.  This  matual  interaction  between  the 
sources  produced  intermodulation  of  the  two  signals  in  the  drivers  them- 
selves. The  weak  signal  drivers  were  particularly  affected  by  the  pump 
signal.  The  branch  of  the  yoke  which  held  the  small  signal  driver  was 
filled  with  size  00  steel  wool  to  change  the  acoustic  impedance  of  the 
branch.  The  steel  wool  decreased,  by  approximately  5 dB,  the  amplitude  of 
the  received  pump  signal  by  the  small  signal  driver. 

3-  The  Receive  System 

The  following  list  of  equipment  comprised  the  receive  system: 

a.  Bruel  and  Kjaer  4-136  l/4  in.  condenser  microphones, 

b.  Bruel  and  Kjaer  2619  field  effect  transistor  (FET) 
preamplifiers, 

c.  Bruel  and  Kjaer  2803  dual  channel  microphone  power  supply, 

d.  Bruel  and  Kjaer  4220  microphone  calibrator, 

e.  Bruel  and  Kjaer  A00029  30  m microphone  extension  cables, 

f.  Hewlett  Packard  3580A  spectrum  analyzer, 

g.  Hewlett  Packard  400  EL  FET  voltmeter, 

h.  Tektronix  545B  dual  channel  oscilloscope,  and 

i.  Honeywell  500  2X-Y  recorder. 
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Figure  4.5  is  a diagram  of  the  receive  system.  A microphone  was  always 
in  the  monitor  port  (the  microphone  port  in  the  driver  mount)  to  continu- 
ously monitor  the  output  of  the  drivers.  To  observe  the  propagation  of  a 
wave  as  a function  of  distance,  a second  microphone  was  alternately  used 
in  the  seven  remaining  microphone  ports.  The  second  microphone  is  called 
the  downstream  microphone.  At  the  beginning  and  end  of  an  experiment, 
both  microphones  were  calibrated  with  a B&K  4-220  pistonphone. 

B.  Test  and  Evaluation 

1.  The  Transmit  System 

A useful  way  to  evaluate  the  transmit  system  is  the  "black  box 
approach."  In  the  black  box  approach,  once  a desired  waveform  has  been 
observed  at  the  monitor  microphone  location,  the  manner  in  which  the  wave- 
form is  produced  is  of  no  consequence.  The  propagation  distances  are  not 
measured  with  respect  to  the  monitor  microphone  position,  but  the  monitor 
microphone  port  is  the  closest  observation  point  of  the  source  waveform. 
Whatever  is  observed  downstream  from  the  monitor  microphone  is  due  to  the 
waveform  observed  there,  or  some  phenomenon  downstream,  and  is  not  due  to 
the  transmit  system.  Valid,  reliable  experimental  results  depend  on  the 
accuracy  of  the  feceive  system. 

Getting  a desired  waveform  at  the  monitor  microphone  location  is 
not  a trivial  problem.  Once  a desired  waveform  is  produced,  it  also  has 
to  be  stable  and  reproducible.  A desirable  waveform  is  an  acceptable 
approximation  of  a pure  tone.  If  the  second  harmonic  of  a tone  is  at 
least  20  dB  smaller  in  -amplitude  than  the  fundamental,  the  waveform  of  the 
tone  is  deemed  satisfactory.  Before  a wave  propagates  to  the  monitor 
microphone  location,  it  experiences  distortion  from  three  sources: 
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(l)  distortion  in  the  transmitting  electronics,  (2)  distortion  in  the 
drivers,  and  (3)  finite  amplitude  distortion  with  propagation.  For  these 
reasons  pure  tones  are  very  difficult  to  produce  at  the  levels  necessary 
for  finite  amplitude  research.  For  the  experiments  reported  in  Chapter  V, 
in  all  cases,  the  second  harmonic  was  more  than  35  <TB  down  from  the  funda- 

4 

mental  in  the  electrical  signals  to  the  drivers.  The  requirement  that  the 
second  harmonic  be  at  least  20  dB  down  from  the  fundamental  at  the  monitor 
microphone  location  was  met  in  all  of  the  experimental  runs  except  one. 

In  this  run  the  pump  had  an  amplitude  of  1^9  dB  and  a frequency  of  3-5  kHz 
and  the  second  harmonic  at  the  monitor  was  only  10  dB  down  from  the  funda- 
mental. The  JBL  375-H  driver  with  an  aluminum  diaphragm  was  found  to  be 
prone  to  strong  distortion  at  3*5  kHz. 

The  amplitudes  required  from  the  drivers  (especially  the 
JBL  375-H  driver  with  an  aluminum  diaphragm)  in  producing  the  desired 
acoustic  signals  often  led  to  fatigue  and  eventual  breakdown  of  the 
drivers.  Evidence  that  a driver  was  experiencing  an  input  overload  was 
the  inconstant  output  amplitude  from  the  driver  and  the  abnormally  high 
second  harmonic.  A necessary  safeguard  to  protect  a driver  was  a quick 
blow  fuse  on  the  output  of  the  amplifier  which  just  barely  permitted 
operation  at  the  desired  energy  requirement.  Another  precaution  was  to 
use  an  input  grounding  switch  on  the  input  to  the  amplifiers.  As  an 
experiment  progressed,  the  downstream  microphone  was  moved  to  various 
locations  to  obtain  waveforms  as  a function  of  propagation  distance.  The 
acoustic  signal  in  the  pipe  was  turned  off  when  the  microphone  was  moved. 
The  signal  to  the  drivers  was  quieted  by  grounding  the  input  to  the 
amplifiers,  and  therefore  the  simplifiers  remained  on  during  the  numerous 
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microphone  location  changes.  Before  the  input  was  either  grounded  or 
turned  on,  the  amplitude  of  the  input  signal  was  decreased  by  at  least 
20  dB  so  that  any  transient  caused  by  throwing  the  grounding  switch  would 
not  be  severe;  large  transients  may  cause  the  voice  coil  of  a driver  to 
exceed  its  maximum  allowable  outward  excursion  and  may  possibly  damage  the 
driver.  The  JBL  375-H  driver  failed  in  two  ways.  A fusible  link  to  the 
voice  coil  melted  (resulting  in  no  response  from  the  driver)  or  the  voice 
coil  assembly  overheated  and  broke  down  (resulting  in  a burning  phenolic 
odor  from  the  driver  mouth) . 


2.  Plane  Wave  Pipe 

In  Chapter  II  it  was  shown  that  the  major  contributor  to 
attenuation  in  the  pipe  was  attenuation  due  to  boundary  layer  effects 
(Eq.  2.37).  Under  normal  experimental  conditions  (temperature  23°C, 
relative  humidity  and  atmospheric  pressure  752  mm  Hg)  the  theoretical 

value  of  is  3.60  X 10-1+  Np/m.  Schaffer^  performed  an  experiment  in 
the  plane  wave  tube  to  measure  0^;  Fig.  4.6  is  a plot  of  the  results.  The 
experimental  value  of  was  found  to  be  3*37  X 10  ^ yJT  - O.OO558  Np/m. 

In  the  theoretical  work  performed  in  this  thesis,  the  experimental  value 
of  0^  is  used. 

The  absorbing  termination  in  the  last  section  of  pipe  was 

critical  in  satisfying  the  progressive  wave  assumption,  especially  for  the 

* 28 
pump  because  of  its  large  amplitude.  A paper  by  Burns  served  as  the 

basis  for  an  experimental  study  of  various  materials  for,  and  shapes  of, 

end  terminations.  A tone  burst  scheme  was  used  to  evaluate  different 

materials  and  shapes  for  the  termination.  Figure  4.7  is  a plot  of  the 

data  for  the  best  all-purpose  termination  found.  The  plot  is  for  a 


FIGURE  4.6 

COMPARISON  OF  KIRCHHOFF  ATTENUATION  THEORY  WITH  OBSERVED  DATA 


(Taken  from  Schaffer^) 


<9  DEMERIT 
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tangent  taper  termination  made  of  fiberglass  insulation  material;  the 
termination  consisted  of  a meter  long  tangent  taper  section  followed  by  a 
solid  cylindrical  section  (see  the  insert  in  the  figure).  Frequency  is 
plotted  along  the  abscissa,  and  the  demerit  or  reflection  coefficient  is 
plotted  along  the  ordinate  in  the  plot.  Demerit  is  defined  as  the  ratio 
of  the  amplitude  of  the  reflected  signal  to  that  of  the  incident  signal; 
the  corresponding  decrease  in  decibels  of  the  reflected  s_0  al  is  ritten 
next  to  the  demerit  values.  A demerit  of  0.01  was  the  s~'gn  goal  (the 
dashed  lin<=  in  the  figure);  the  fiberglass  tangent  taper  provided  a 
demerit  of  at  least  0.01  for  frequencies  above  1250  Hz.  In  all  the  experi- 
mental runs  the  frequency  of  the  pump  was  greater  than  1250  Hz. 

This  thesis  concerns  the  interaction  of  plane  waves;  under  what 
conditions  is  the  propagation  of  acoustic  signals  in  a pipe  planar? 

Planar  propagation  in  pipes  and  ducts  is  certain  only  for  frequencies 
below  a cutoff  frequency.  The  cutoff  frequency  for  the  various  propaga- 
tion modes  in  a circular  tube  of  radius  r is  expressed  by 

D c 

f _ gs  ° 

mn  2itr  ’ 

where  is  the  nth  root  of  the  following  equation : 

J '(D  ) = 0 

The  lowest  cutoff  frequency  is  the  one  for  m=l,  n=0.  This  frequency  is 
59^0  Hz  for  the  tube  used  in  the  experiments  reported  here.  A signal  of 
frequency  greater  than  3960  Hz  will  not  necessarily  propagate  in  a non- 
planar  mode,  and  it  was  generally  found  that  experiments  with  primary 
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frequencies  above  the  lowest  cutoff  were  run  apparently  without 
stimulating  nonplanar  modes.  If  a fundamental  undergoes  planar  propaga- 
tion, the  higher  harmonics  of  the  fundamental  that  have  frequencies 
greater  than  the  lowest  cutoff  frequency  will  nevertheless  also  tend  to 
have  planar  propagation,  because  the  higher  harmonics  are  produced  by  the 
propagation  of  the  fundamental.  The  exception  to  this  condition  is  pro- 
vided by  irregularities  in  the  pipe,  such  as  a poor  flange  connection  and 
a protruding  microphone  or  plug,  which  create  a boundary  condition 
different  from  the  original  one.  These  new  boundary  conditions  require 
different  modes  of  propagation  of  the  signal  to  satisfy  them.  Thus,  the 
pipe  irregularities  scatter  energy  out  of  the  plane  wave  mode  into  higher 
order  modes.  If  the  irregularities  are  small,  the  amount  of  energy  lost 
from  the  plane  wave  mode  is  small. 

The  actual  propagation  distances  corresponding  to  the  microphone 
locations  are  listed  in  Table  4.1. 

TABLE  4.1 

DISTANCES  FROM  SOURCES  TO  MICROPHONE 
MEASURING  LOCATIONS  IN  METERS 


MIKE 

STATION 

MONITOR 

n 

2 

3 

D 

6 

7 

IHQH 

4.2 

m 

11.6 

15.3 

19.0 

22.8 

26.5 

DOUBLE  YOKE 

0.1 

3.8 

7.6 

11.3 

15.0 

18.7 

22.4 

26.1 

The  propagation  distances  are  measured  with  respect  to  the  finite  amplitude 
producing  drivers.  When  the  single  yoke  was  used,  the  propagation 
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distances  were  measured  with  respect  to  the  face  of  the  driver  producing 
the  pump  signal.  The  propagation  distances  for  the  three-driver  pump 
source  (the  double  yoke)  were  measured  with  respect  to  the  point  of  the 
yoke  where  the  three  individual  pump  signals  merged.  The  source  ampli- 
tudes for  the  two  signals  (pump  and  weak  signal)  were  obtained  through 
interpolation  (by  the  use  of  the  Fubini  solution  or  small-signal  theory) 
of  their  amplitudes  observed  at  the  monitor  microphone  location. 

J.  The  Receive  System 

The  receive  system,  treated  as  a single  unit,  was  tested  for 
distortion  of  the  received  signals  and  needed  to  be  calibrated.  The 
distortion  test  proceeded  as  follows.  The  microphone  cartridge  was 
removed  from  the  preamplifier  of  the  B&K  l/^  in.  microphone.  Two  elec- 
trical signals  of  equivalent  frequency  and  amplitude  as  produced  in  an 
experimental  run  were  introduced  to  the  preamplifier  through  an  electronic 
signal  summer.  The  spectrum  of  the  signal  produced  was  viewed  on  the 
spectrum  analyzer  for  intermodulation  distortion  components  located  at 
frequencies  where  sidebands  of  the  acoustical  interaction  would  be.  The 
results  for  the  equivalent  input  of  a signal  of  4. 5 kHz  at  a SPL  of  150  dB 
and  a signal  of  1 kHz  at  a SPL  of  110  dB  were  that  no  intermodulation 
distortion  component  had  an  amplitude  larger  than  TO  dB.  Since  the  ampli- 
tude of  the  sidebands  recorded  during  an  experiment  were  usually  around 
100  dB,  they  were  therefore  not  produced  in  the  receive  system.  This 
distortion  test  checked  the  distortion  of  the  complete  receive  system 
except  for  the  B&K  microphone  cartridge.  It  was  believed  that  the  dis- 
tortion caused  by  the  cartridge  would  be  smaller  than  the  distortion 
caused  by  the  rest  of  the  system.  The  calibration  of  the  system  was  begun 
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by  calibrating  the  spectrum  analyzer,  according  to  its  manual,  for 
frequency  and  amplitude  response.  The  second  part  of  the  procedure  was 
to  calibrate  the  entire  system  with  the  B&K  pistonphone.  The  receive 
system  was  calibrated  before  and  after  every  experimental  run. 


CHAPTER  V 


COMPARISON  OF  EXPERIMENTAL  RESULTS  AND  THEORETICAL  PREDICTIONS 


A representative  sample  of  the  experimental  results  for  Case  II 

are  presented  here  and  compared  to  the  theory  developed  in  Chapter  II. 

For  the  weak  signal  the  experiments  covered  a range  of  source  SPL  from 

104  to  121  dB  and  a frequency  range  from  0.6  to  1 kHz.  The  corresponding 

ranges  for  the  pump  were  125  to  158  dB  and  1.5  to  6.6  kHz,  respectively. 

The  measurement  distances  ranged  from  3-9  m to  27  m.  The  largest  value  of 

Op  for  which  a measurement  was  made  was  30-  A total  of  25  independent 

experiments  were  run.  Because  of  the  similarity  and  conformity  of  the 

results,  only  three  of  the  experimental  runs  are  presented  here.  No 

experimental  waveforms  are  presented.  The  relative  amplitudes  of  the  weak 

signal  and  the  pump  were  such  that  the  width  of  the  oscillogram  trace 

overshadowed  the  modulation  of  the  pump  by  the  weak  signal  (the  modulation 

may  be  seen  in  the  theoretical  predictions,  Figs.  1.2  and  3-5) • 

An  example  spectrum  of  an  experimental  waveform  is  shown  in 

Fig.  5.1.  Spectra  for  u ',  u ',  and  u are  shown  in  plots  (a),  (b),  and 

w p 

(c),  respectively.  The  ordinates  are  SPL  and  the  abscissas  are  frequency. 
The  values  of  the  indices  above  a component  identify  its  relative  fre- 
quency (or  its  u representation).  The  notation  is  that  presented  in 
m,n 

Chapter  II.  Some  components  are  the  result  of  degeneracy,  i.e.,  the 
addition  of  two  or  more  components  whose  frequencies  are  the  same. 
Distortion  of  the  weak  signal  may  be  seen  in  plot  (a)  (the  second  harmonic 
of  the  weak  signal  is  26  dB  down  from  the  weak  signal).  The  "spurious” 


SPL  - dB 


95 


weak  signal  and  sidebands  in  plot  (b),  which  should  be  the  spectrum  of  the 
pump  signal  by  itself,  are  caused  by  the  presence  of  the  weak  signal  in 
the  signal  to  the  pump  drivers.  The  weak  signal  leaked  into  the  pump 
transmit  system  even  though  the  input  to  the  weak  signal  amplifier  was 
grounded.  Leakage  of  the  pump  signal  was  avoided  in  plot  (a)  by  grounding 
the  input  of  the  pump  amplifier  and  disconnecting  the  pump  signal  from  the 
grounding  switch.  The  sidebands  are  caused  by  the  interaction  of  the  weak 
signal  and  its  harmonics  with  the  pump  and  its  harmonics.  The  -1,1  com- 
ponent is  the  first  lower  sideband  and  the  1,1  component  is  the  first 
upper  sideband.  The  dashed  lines  above  the  1,0  and  o,n  components  in 
plot  (c)  represent  the  SPL  of  the  corresponding  components  in  the  other 
two  plots.  The  difference  between  the  dashed  line  and  the  amplitude  of 
the  components  in  plot  (c)  is  1 dB,  1 dB,  1.5  dB,  2.2  dB,  and  3*2  dB  for 
the  components  1,0  and.  o,n  (n=l,2,3^)j  respectively.  Little  decrease  of 
the  weak  signal  is  observed,  but  the  pump  and  its  harmonics  are  definitely 
reduced.  The  higher  the  harmonic  of  the  pump,  the  larger  the  difference 
in  SPL  of  the  harmonic  in  plots  (c)  and  (b)j  this  trend  was  predicted  in 
Chapter  II  (Eq.  2.26).  The  difference  between  the  two  SPLs  is  a measure 
of  the  amount  of  energy  redistributed  from  the  harmonics  to  the  sideband 
formation. 

Figure  5.2  shows  a comparison  of  predictions  and  measurements  for 
an  experiment  in  which  the  source  conditions  were  as  follows: 

SPL  = 1^0  dB,  SPL  = 115  dB,  f = 1-5  kHz,  and  f =1  kHz. 

p 7 w 7 p w 

The  changes  in  SPL  of  the  weak  signal  (u^  Q) , the  first  lower  sideband 
(u  ,),  the  pump  (un  .),  and  the  first  upper  sideband  (u.  ,)  as  a 

"1*1.  O • J.  -*- » J- 
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(b)  FIRST  LOWER  SIDEBAND  (u.,,  , ) 


0 1 2 


(c)  PUMP  (u0>  , ) 


FIGURE  5.2 


COMPARISON  OF  EXPERIMENTAL  DATA  WITH  THEORY 
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function  of  ct^  are  shown  in  plots  (a),  (b),  (c),  and  (d),  respectively. 

In  plots  (a)  and  (c)  the  dashed  lines  represent  the  decrease  of  the 
respective  components  based  on  small-signal  absorption.  The  difference 
between  the  small-signal  predictions  and  the  weak  shock  theory  predictions 
for  the  weak  signal  represents  the  suppression  of  the  weak  signal.  In  the 
case  of  the  pump  [plot  (c)]  the  difference  primarily  represents  the  excess 
attenuation  caused  by  the  finite  amplitude  distortion  of  the  pump  itself. 
The  particular  experimental  run  shown  in  Fig.  5.2  is  a little  unusual 
because  the  first  lower  sideband  is  to  the  left  of  the  weak  signal  in  the 
frequency  domain.  The  two  primaries  were  chosen  close  in  frequency  so 
that  the  value  of  ft  would  be  near  one.  In  Chapter  II,  in  the  analysis  of 
the  Case  II  preshock  solutions,  it  was  observed  that  the  greatest  predicted 
suppression  of  the  weak  signal  occurs  when  ft  has  a value  of  one.  The 
agreement  between  theory  and  experiment  is  good;  the  standard  deviation  of 
the  data  about  the  theoretical  curve  has  a value  of  0.7  dB.  Only  0.8  dB 
of  suppression  of  the  weak  signal  was  achieved.  The  first  upper  sideband 
grows  more  rapidly  and  is  larger  than  the  first  lower  sideband.  Such  an 
asymmetry  agrees  with  the  theoretical  predictions  in  Chapter  II. 

Figure  5.3  shows  the  comparison  of  predictions  and  measurements 
for  an  experiment  with  a smaller  value  of  ft.  The  source  conditions  were 
as  follows: 

SPL  = ikb  dB,  SPL  = 113  dB,  f = 2.5  kHz,  and  f =1  kHz. 

The  weak  signal  does  not  undergo  much  suppression  and,  in  fact,  the 
suppression  seems  to  taper  off  at  approximately  1 dB,  around  a^=2.  From 
that  point  on,  the  predicted  decrease  in  the  weak  signal  parallels  the 
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FIGURE  5.3 

COMPARISON  OF  EXPERIMENTAL  DATA  WITH  THEORY 

SOURCE  SPLW  = 1 13  dB,  fw  = 1.0  l<H« 
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linear  theory  curve.  The  agreement  between  theory  and  experiment  is  in 
general  better  than  for  the  previous  experiment;  the  standard  deviation  of 
the  data  about  the  theoretical  curve  is  equal  to  0.5  dB.  Again  the  first 
upper  sideband  is  seen  to  grow  more  rapidly  and  to  be  larger  in  amplitude 
than  the  first  lower  sideband.  At  the  larger  values  of  a , however,  the 
greater  attenuation  of  the  higher  frequency  sideband  has  its  effect  and 
the  first  lower  sideband  becomes  the  larger  one.  The  squiggle  in  the  weak 
shock  theory  curves  is  partially  due  to  the  polynomial  fit  performed  in 
program  GRAPH. 

The  majority  of  the  experimental  results  and  theoretical 
predictions  were  similar  to  Figs.  5*2  and  5*3;  an  exception  is  shown  in 
Fig.  5*4.  In  this  example  (which  shows  the  smallest  value  of  Cl)  the 
source  conditions  were: 

SPL  = 145  dB,  SPL  = ll4  dB,  f =4  kHz,  and  f =1  kHz. 
p 7 w p w 


What  sets  this  particular  experimental  run  apart  is  the  fact  that  the  weak 
signal  undergoes  amplification  (relative  to  the  small-signal  prediction) 


rather  than  suppression.  The  enhancement  of  the  weak  signal  (approximately 


1 dB)  appears  to  occur  between  the  values  a =1  and  a =4.  Beyond  the  point 

P P 

0^=4,  the  weak  signal  decays  according  to  linear  (small- signal)  theory. 


The  agreement  between  theory  and  experiment  is  good;  the  standard  devia- 
tion of  the  data  about  the  theoretical  curve  is  equal  to  0.7  dB.  Again  the 
upper  sideband  is  more  prominent  at  first;  the  lower  sideband  is  more 
prominent  later.  An  extensive  experimental  and  theoretical  investigation 
was  undertaken  to  study  the  enhancement  of  the  weak  signal,  but  the 
results  of  the  study  were  not  conclusive.  The  theoretical  work  showed  that, 
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whenever  enhancement  occurred,  the  pump  frequency  was  a harmonic  of  the 
weak  signal;  the  investigation  also  showed  that  this  frequency  condition 
was  necessary  but  not  sufficient.  The  predicted  enhancement  was  amplitude 
dependent  as  well  as  frequency  dependent.  The  experimental  work  was  less 
edifying.  No  other  experimental  run  demonstrated  the  enhancement  of  the 
weak  signal  to  the  degree  seen  in  Pig.  5 Figure  is  presented  to 
show  a result  which  was  observed  both  experimentally  and  theoretically 
but  which  is  not  completely  understood. 


CHAPTER  VI 


CONCLUSIONS  AND  FINAL  COMMENT 

The  purpose  of  this  research  has  been  to  6tudy  the  suppression 
of  sound  by  sound  of  a higher  frequency.  More  specifically,  the  planar 
propagation  in  air  of  a finite  amplitude  tone  (the  pump)  and  a smaller 
amplitude  tone  (the  weak  signal)  of  lower  frequency  has  been  investigated 
both  theoretically  and  experimentally.  A brief  summary  of  the  work  com- 
pleted and  the  conclusions  drawn  from  those  results  follows. 

Suppression  is  caused  by  a modulation  of  the  weak  signal  by  the 
pump.  The  acoustic  energy  of  the  weak  signal  is  not  lost  (absorbed)  but 
is  rather  redistributed  to  other  parts  of  the  spectrum.  The  amount  of 
suppression  depends  qualitatively  on  whether  the  pump  frequency  is  lower 
than  the  weak  signal  frequency  (Case  I)  or  vice  versa  (Case  II).  In 
Chapter  II,  Case  I and  Case  II  interactions  were  reviewed  and  compared. 

A general  notation  was  introduced  so  that  the  various  formulas  could  be 
used  for  either  Case  I or  Case  II.  The  exact  solution  valid  before  shock 
formation  was  reviewed,  both  for  the  time  domain  (Eamshaw  solution)  and 
for  the  frequency  domain  (Fenlon  solution).  A convenient  way  of  writing 
Fenlon's  solution  was  given.  Various  approximate  solutions  were  formed 
from  the  exact  solutions  and  compared  with  each  other;  analysis  of  these 
solutions  led  to  many  conclusions.  Suppression  before  shock  formation  in 
Case  I may  easily  exceed  40  dB,  but  for  Case  II  suppression  may  never 
exceed  2.4  dB  and  in  most  practical  cases  is  much  less.  The  modulation 
process  that  causes  suppression  also  leads  to  the  formation  of  sidebands 

Lk 
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about  the  high  frequency  signal.  The  primary  energy  source  for  the 
sideband  formation  in  both  Case  I and  Case  II  was  observed  to  be  the 
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higher  frequency  signal  and  its  harmonics,  if  any.  In  Case  II,  it  was 
observed  that  the  higher  the  harmonic  of  the  pump,  the  greater  the 
harmonic's  amplitude  is  decreased  in  the  interaction.  The  first  upper 
sideband  grows  more  rapidly  (at  least  at  first)  and  is  larger  than  the 
first  lower  sideband. 

Theoretical  predictions  valid  beyond  the  shock  formation 
distance  were  obtained.  Included  in  these  predictions  are  the  effects  of 
small-signal  attenuation  and  dispersion.  Weak  shock  theory,  cast  in 
difference  equation  form,  was  used  to  obtain  the  predictions.  In 


Chapter  III  the  difference  equations  were  implemented  in  a computer  algo- 
rithm. The  computer  program,  referred  to  as  program  C500,  incorporates 
attenuation  and  dispersion. 

The  agreement  between  the  experimental  results  and  the 


theoretical  predictions  shown  in  Chapter  V was  excellent  and,  in  general, 
supported  the  conclusions  drawn  from  the  pre shock  solutions.  The  experi- 
ments, which  were  for  Case  II  only,  showed  little  suppression  of  the  weak 
signal;  the  maximum  suppression  observed  was  about  1 dB.  The  suppression 


of  the  weak  signal  occurred  at  relatively  small  values  of  (a  is  the 
shock  formation  distance  of  the  pump);  for  values  of  cr^  approximately 
equal  to  3>  the  suppression  process  halted  and  the  weak  signal  propagated 
thereafter  according  to  linear  theory.  The  first  upper  sidebands  were 
observed  to  be  more  prominent  than  the  first  lower  sidebands  for  small 


values  of  a but,  because  of  absorption,  for  larger  values  of  the  lower 
sidebands  were  more  prominent.  In  some  cases  the  weak  signal  became 
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enhanced  (relatively)  rather  than  suppressed.  Enhancement  was  observed 
both  experimentally  and  theoretically;  it  appeared  to  occur  when  certain 
frequency  and  amplitude  conditions  were  met  by  the  primaries.  In  particu- 
lar, the  pump  and  weak  signal  had  to  be  harmonically  related.  With  the 
primaries  so  related,  degeneracy  occurs.  Thus  the  enhancement  may  have 
been  caused  by  the  presence  of  sidebands  equal  in  frequency  to  that  of  the 
weak  signal. 

The  primary  task  undertaken,  the  investigation  of  Case  II,  was 
accomplished.  The  conclusion  is  that  little  suppression  of  a weak  signal 
by  a higher  frequency  pump  is  possible,  before  or  after  shock  formation. 

It  has  been  firmly  established  here  that  if  large  suppression  of  a tone 
is  desired,  the  suppression  may  only  be  obtained  through  a Case  I inter- 
action, even  though  the  Case  II  interaction  is  more  attractive  in  respect 
to  harmonic  and  sideband  "pollution."  The  computer  programs  developed  to 
implement  the  theoretical  investigation  of  Case  II  are  applicable  for 
other  purposes.  The  computer  programs  may  easily  be  adapted  to  other 
types  of  propagation  (spherical  or  cylindrical)  and  to  other  forms  of 
attenuation  (atmospheric  or  underwater). 

The  present  work  completed,  what  is  the  next  step?  A logical 
extension  of  the  work  presented  here  is  the  substitution  of  noise,  either 
narrowband  or  wideband,  for  one  or  both  of  the  two  acoustic  signals.  The 
work  with  tones  would  seem  to  indicate  that  the  most  efficient  way  to 
suppress  a band  of  noise  would  be  a Case  I interaction  with  a lower  fre- 
quency pump.  The  theoretical  work  in  this  case  could  then  be  restricted 
to  shock- free  regions.  The  preshock  solutions  are,  in  general,  easier  to 
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apply  than  the  postshock  solutions.  In  any  case,  however,  the  applica- 
tion to  suppression  of  random  noise  is  not  a trivial  theoretical  problem. 

g 

Some  work  involving  noise  has  already  been  done.  Pestorius  has  studied 
the  propagation  of  wideband  finite  amplitude  noise  (no  tones  present). 

I 

Some  theoretical  work  on  the  interaction  of  noise  with  a tone  has  been 
reported  by  Rudenko,  Soluyan,  and  Khokhlov.^0  Work  is  presently  being 
done  at  ARL  on  the  interaction  of  a finite  amplitude  tone  and  a bandpass 
of  noise  (Ref.  26) . 

| 


APPENDIX  A 

LISTINGS  OF  PROGRAM  C500  AND 
SUBROUTINES  USED  BY  PROGRAM  C500 


* 
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PROGRAM  CSOO 

COMMON  T <*?•«!  .1114260)  .»  T I I BOO) 
DIMENSION  |IU200>  .LU<M2l 
EQUIVALENCE  iu. Ill 

REAL  MAI]. MAI 2 


fma*b?oooo. 

MRHF  12)  A.F|.R.r?ttMrTA.OISIANCCtUlSlMCtMRt).FMAA,MO 
LPPlHTSaNHCNO 


C THIS  IS  PROGRAM  CSOO  MMICM  PROPAGATES  AM  iMRut  ■ A Vt I ORM  BT  INC**» 
C MENTAL  STEMS  TO  A OCSIRCO  UISTanCE  US|M«  bIak  SHOC*  »«tOR*.  a" 
C fMUATION  AMO  OlSPfMSIOM  ARE  APPLIED.  » 0»  A Of  T A J lEO  EAPlANATJOM 
C SEC  AML-TR-TT-B  by  RILL  MlLLSHiRt.  The  FOtlOvING  IS  A LIST  OF 
C THE  IMPORTANT  PROGRAM  MARAMEIEm$. 

c a-spl  of  small  signal  re  .00002  ht/m 

c fi-frcoucncy  or  spall  siomal 

C B-SPL  or  HIGH  MCOOCNCV  pump 

c fi-  rRfouCNCt  or  pump 

C DISTANCE-  OESiREO  PROPAGATION  ulSIANCE 

C 01 S INC-  INCREMENTAL  PROPAGATION  UlSTANLE 

C NC ip?-  TmC  mo.  or  CLEMENTS  IN  U A no  T ARRA1CS  ■ ?«*NC1P2 

c rpAi-  rarouCNCT  or  last  fouriem  cocff.  calculated  bt  fft 

C SCAL-  USCO  TO  SCALE  H ARRAY  BEFORE  FFT 
C U-  PARTICLE  VELOCITY  ARRAY 
C T.  CORRESPONDING  TIME  ARRAY 
C vSCALaSCAL'  USED  FOR  INI T I AL I ZAT ION 
C FOUMmF  MAI • USED  FO*  | N | T I AL 1 2 A I I 06- 

C~  MATPIS-NO  or  PROP.  STEPS  BCTaCEN  ATTEN  AMO  OISPER  CORRECTIONS  * 

C NATOIS*NOISTOMS 

C MATOISRNATOIS.  USED  FOR  initialization 

c Na?«*NElP?t  The  NO  or  ELEMENTS  IN  U MV  t ARR AIMS  INI T 1ALL V 
C NPOINTS-  NO  Or  ELEMENTS  IN  V MO  1 AFTER  CALL  TO  SUB.  RAVCRAOP 
C NPOINTS-  USCO  TO  INITIAL.  NPOINTS  mhEN  ATT-OISP.  NOT  APPlIEO 
C LPOINTS-  USED  TO  INITIAL.  NPOINTS  MHEN  ATT-OIS  ARC  APPLIED 
C NMMN/? 

C NM| a 1 «N/?  . IS  The  no.  or  FOURIEM  COEFFS.  calculated  bt  fft 
c NRO-NO  or  MUMP  RAVE  lengths  in  initial  U ARRAY 
C OELT-  ORIGINAL  RAVEFO0M  I NCRCMEN I 
C MOO-  NO  OF  ZEROES  In  RIGHT  KRO  BUFFER 
C NOOL-  NO  OF  ZEROES  IN  LEFT  ZERO  RUFFE* 

C NGOaNOOL • 1 * FIRST  NON-ZERO  ELE . IN  U ARRAT 
C NENOaNPO I N I S-NOOL • LAST  NON-ZEMO  El£.  IN  U array 
C NROOaNCNO • I • FIRST  ZC»0  ElE.  IN  R|6mT  2E«0  BUFFER 
C NflEN0aNCND*M00*  LAST  ZERO  ELE.  In  u ARmay  In  RIOmt  ZERO  SuFFCR 
C NOISTORS-  NO  OF  InCRE • PROP.  S1EPS  PERFORNt  by  ONE  call  TO 
C SURROUTINE  MAVEPROP 

C NO-NO  nr  CALLS  TO  RAVtPROP  TO  REACH  UCSI*EU  PROP.  OISTANCE 
C 

•aaaaaaaaaaaaaaaaaaaaaaaaaaaaaa •••#•• 


C a aaa a* •••##*«*•*••*•• •••••*•• 

•EAO  ZSA.A.Fl.i.FJ.THf TA 
2S0  F0RMATI5FB.2) 

READ  260. OISTANCE. DISTNC.NfiPZ 
2*0  FORMAT (2FR.2* IS) 

•EAO  2T0.Fmaa.SCAL 
2*0  format Ljr io. el 
PlaJ.lalS* 


C AMR  OF  TONES  CAST  INTO  PARMCLt  VtL.  FmQm  >pl  Rt  .00002  NT/m 2 


AAa( .0000000*B2>*lO.#aiA/20.» 
68a t . 0000000*02 1 *10.** <0/20.1 

■SCALaSCAL 

F pi  )MaF  M A I 

NATO  I S»2 

NOISTOMSa2 

MATU I 5aNA  T o I S 

NPO|NTSa2.0a*'NFIP2) 

N.NPD I NTS 

LPOlNTSaNPniNTS 

NMaN/2 

Nh I • 1 N/2 • • 1 

NPOar2* IN- l ) / IfaFRAII 

PRINT  a«,NP0 

A?  FORMAT!  • NPO  a a.IlO'. 

OE  L T aNPO / <F2* (NR0INTS-) I » 

NOOaAS 

NOOLalS 

NGOaNOOL* I 

NENOaNPO I NT S»NOOL 

NROOaNPOfNTS* I *NOCl 
NBFNPaNPOINTS.NOOL  *NOn 


r DL  101  CREATES  LEFI  ZERO  RUFFE M 


00  101  Jal.NOOL 

U< J) aO. 

T<  J) ajapEtT -DEL  T 


00  601*  Ilal.NC 

01STakCI*CIST an/E *0 1 5 1 NC*NO 1ST  ORS 

Call  fftcm(i.u.aa*o.at» 

Call  fetcmiijs.t.Iavo.aT) 

CALL  01 SKPA I T 
SCALaaSCAL 

S I G*  J . 200**80*2 . aP I *F  2*0 1 STANCE / <J*0.*JA0. I 
NPOlNTSauPPlNTS 


C SUB  MAVEPROP  PROP Ab A TE S THE  MAVCFORR  VIA  FEAR  $MOC*  THEORY 

c 



CALL  *A VC  PROP (U.T.O! ST ANCE.01SINC.NP0INTS. NOISTORS) 

NPrlNTS*NROINTS 
Call  store < l .u.baoo.a?) 

Call  STORE l 133.T.0*0«.*T) 

Call  oisrmait 

C..........M........a.....aaa...ia....a....a.........a....i.ta.... 

c 

C SUB  TINE  frc  REMOVES  The  ZERO  buffers  and  peRFOhnS  a fft 
C 

C •••••••••••••*.#•••••**••••••••*•**•******•••••*••**••***• .a. aaaa* 

CALL  TINEFRE IU.T.NP0INTS.N.NEIP2.II.SCAL.NM. I TOT. IMAA.LU.I I .SIG. 
CNAfilF  .OT.NRENr!) 

C..................a.a.**«...M....a.aM«..M....aaa....**.a.t.i... 

C 

C SUP  HANGCN2  APPLIES  A HANNING  UA 1 A ■ 1 NUO*  ’0  THt  SPECTRUM  OF  W 
C INTERMEDIATE  MAVEFORM 


Call  mangEnein.td 
222?  SCAL*SCALa>.vaaFLCAT ' IT0T) 


C OL  20  SCALES  TmC  SMECTRum  COmPUNENTS  Tu  ABSOLUTE  VALUES 


DO  20  L*|.M 
jbn.J-L 

U I J • a 1 1 I J I / <N*SC*L > 
20  CONTINUE 
ShORbl '2> 


"L  31S  CALCULATES  The  FOURIER  COEFFICIENTS 


00  315  |b2.NN 
J*?«t!-l>*l 
JJaJ* I 

AmP«SORT IU< J)*U< J) *U< JJ)*U< JJ) ) 

IF (Amp. CO. 0. ) GO  TO  mb o 

U I I I a?0.**LOb!0  I ( A l 5*  AMP  I / I < • VOO  V? ' • I 1 .*!«>  » I 
GO  To  69? 

869  U ( I I aO. 

690  CONTINUE 

IF <u<I I .LT.O.Ol  UlliaO.U 
3 is  continue 

IF  IUIJ I .CO.O.)  GO  TO  «01 

U<  I ) a20 . * ALOG 10(AlS*ABS<U<Ii I / < 1 . A U* . U0U02 ) » 
GO  TO  692 
69)  U(1>*0. 

89?  Continue 

IF <u<l ) .LT.O.O)  U ( 1 i *0 . V 
IF (SHOP  .EC. v.  i oc  TO  «93 

UINMI I *26.aALOGlO<A|SaABS  < SHU* ) / I 1 .*!«•. U00C2>  > 
GO  TO  694 
R9J  UINMI I *0 . 

§9*  CONTINUE 

IF (UINMI l .LT.O.O*  U'NMJlao.O 


C OL  666  CALCULATES  l«F  FREQUENCIES  OF  The  Fui'Rlt"  COEFFICIENTS 


DO  886  jal.NMI 
T(J)b(J-D/TNAI 
88P  CONTINUE 
UM*-i 
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ft  i U • 

•o  to  rg 

rt  I l (.Ml  • I I-2«0l /THAI 
FT  I li**n  *u  < i«l  > 

•6  CONTINUE 

s continue 

t\9  CONTINUE 

RR|NT  T01.n.!L» 

to i Tn«***r<  • THf  i sot i is  . 

»P!tr  l?»  SI0tIl*aPT 

LBOtNTSaMBOINTS 


It*  ■ *• *»» 


I»  mTT-OIS  A«r  to  be  APRlUO  OL  MIS  APPLIES  THE"  IF  NOT  CONTROL 
UPTURNS  TO  TmI  BEGINNING  10  UL  M010 


c«« 


C TlHCPRf  REMOVES  ZC«0  BUT  FEUS  ANQ  PERF  OHMS  TP  T ON  MAVCFOftM 


C O L 01S  APPLIES  0 T TtNUA  MON  ANO  DISPERSION  CORRECTIONS 


00  MIS  1 bJ.NH 

Jp>P|I*||a| 

PR(MI1*| l/TNAI 

'0T(*  OHrta  | ,601371  *SORT  IFRE ) - • OOSS  7 7 I *N0 1 ST0NS*0 I S 1NC*MATPI S 
ZZbIiiji 

M I J) aEXP l-ARO) • • I* I J) *C0S I AHGi »lX ( JJ) aSlNIAROi 1 

1 AT  JJ)  al*P  (-ARG)*  • !*«JJ)»COS  lARbi -2Z«SINIAR6>  ) 

■IS  CONTINUE 

IOPaJSC-0 


C««< 


C PERFORMS  BACFbARO  FFT  TO  GE 1 M«CK  INTO  W TI-E  DOMAIN 


lTOTRalPPTOMIIX.NElP2.IOP.ISC.LOi 


C OL  Rift  SCALES  BAVETORN  TO  ABSOLUTE  VALUES  AFTEM  BACKBAHQ  PF f 


00  «)ft  Jal.N 
L»N«1-J 

UIL!»I*IL>/<N«MSCAL**.0«*l!T0T.! |0>M) i 
Rift  CONTINUE 


C OL  BIB  SHIFTS  CLP.  OP  U ANO  T TObARO  TU  hake  ROUn  FUR  LEFT 
C ZERO  RUFFE". 


C ....aa.aa. ....... a. 

00  B2I  lal.NCOL 
UfTlaft. 

T»I>all-ll #0T 
•21  CONTINUE 

LPOlNTSaNBPNO 

HBuFetPOINTS-NOO 


C OL  BIT  CREATES  RIGHT  ZERO  BUFFER 


00  BIT  jaNNQO.NRENO 
UUIbO. 

TlJial  IN)  • I j-NRfiO*  ) l *0T 
•IT  CONTINUE 

CALL  STORE  II .U.RAOO.ATl 

Call  store  iuz.t.maw.aT) 
call  oiskmait 
NaTOI s*nato i s*nato 1 S 

LPOlNTSaN* ] OB 
AO  1 ft  CONTINUE 

RR1NT  l, SI© 

1 format i a tmp  last  sio-a  is  * ••fio.ji 


ENT 


SUBROUTINE  bAvEPROPIU.T.OISTANCE.OISINC.NPOINTS.NUISIORSI 

OlNENSlON  KID. Tin 

COMMON  DUPIB20CI .TSII200) .1  OCSlIFOO) 


r SFC  ARL  TR-TJ-/J  FOM  OCTAIlS  ON  THIS  SUBROUTINE 


BAVEPNOP  is  The  bEak-SmOC*  «*vE  PROmAGAT iun  AlGUR|Thum, 


Call  parameters  — 


U IS  A REAL  ONE  OINtNSIUNAL  AHRAt  UtSCRJBINO  TM|  PART  1 C AL 
BELUClTT  OF  A POINT  on  lHt  BAVE  FOR-  ]N  EITHER  tNGLlS- 
0R  METRIC  UNITS.  IF  MfIRIC  UNI  IS  AME  USEO  Cli.  THE 
SIGNAL  SPEED  OF  SOUNU.  NUS I ME  SUPPLIEO  IN  METRIC  UNITS. 


T IS  A REAL  ONt  OINENSIUNAL  AMRAV  INSCRIBING  THE  TIME  l|N 
SECONDS)  COORDINATE  FOR  A PAT 1CULAR  POINT  ON  THE  a A V£F  0*— . 
The  DOUBLET  (Kill. Till*  Then  DESCRIBES  A bAVEFuhm  POINT 
IN  ThE  PARTICAL  velocitt-tine  PLANE. 


IF  1 1 1 .NE .NATO IS)  on  To  DOlft 
CALL  FrTCH(l.u.RAOw.AT) 
call  fftpmi ij3,t«b*ro,*ti 
CALL  OISKBAIT 
NPO iNTSaMPO I NTS 
BSCALaSCAL 


DISTANCE  -- 

01  STANCE  — A REAL  VARIABLE.  IT  IS  T-t  OISIANCt  OVER  bhICh  The  bavE 
IS  TO  BE  PROPAbATEO. 


OISINC  — A REAL  VARIABLE.  IT  IS  THE  INCHEHENTAL  STEP  SIZE  USEO  TO 
COVER  THE  DISTANCE.  OlSIANCE. 


IROINTS  — AN  INTEGER  VARIABLE.  I I IS  THf  NUMBER  OF  POINTS  US£0  Ti 
OE  SLR  I RE  The  bAvE  PONm  10  Ml  PRURAGATCO. 

C NOISTCRS-  THE  NUmhER  OP  INCREMENTAL  STEPS  PRORAbA TEU  ONE 
C CALL  TO  BAVERROR. 


Call  TIMEFRE IU.T.NR01NTS.N,nE(P2.Ia.SCAL.NH.|TOT. imAA.LU.II.SIG. 
CMAGIF.OT.NRENO) 


PERJOO  a TinPOInTSI 
R ■ 1.205 
COaJAO. 

ETA  a B/C0«a2 
UISF ACTRaETAaOISINC 
00  1 4 0 NOISTuRaj .NOISTOHS 
NS  a nShIM  . D 
TLOB  a TlD-1,0 

CrrCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCLCCCCCCCCCLCCCCCCCCCCCCCCCCCCCCC 


This  loop  tests  fob  Shocks  anc  s«ipts  out  points 

IN  The  OOURLtO  VALUED  RLG10N 


CCCCCCCCCCCCCCCCCCCCCCCCCLCCCCCCCCCLCCLCLCCCLCCCCLCCCLCCCCCCCCCCLL 
00  30  I NOE Aa i .KPOINTS 
NPOlNT  a *Ntr »-N«-lFT 
T.NPOINT)aTllNOE«)-UlINOC«)BOISEACIN 
L • NROI NT ) bu  1 1N0E  * • 

IF  IT (NPOTNT' .mT.TlOb)  GO  10  JO 


TAVG  a 1T1NF  INTI .T iNROlNTin i/2. 0 
DC  10  I a 1 , NS 

IF  ITAvG.LE  • TS  * I ) * GO  10  ill 
CONTINUE 
NS  • NS • I 
TS(NS)  a TAVG 
LCCS'NS)  a NPOlNT 
GO  TO  JO 

NSMIFT  a nSmIFI.NPOTNT-LOCSIII 
LOC  • LOCSUI-I 
T I LOC • 1 » a T(NPOINT) 

U I LOC • I I a UINPOINTl 

TSlD  • IT INROINT) *T ILOC I • /|.0 


NS  ■ I 


f ....... aaaaa.**aa..... a... aaaaa...*. a. a 

00  BIG  lal.h 
NIaN-IM 
NIIaRI.NOOL 
UINII)BU(NT) 

TiMDbT  IND*OT*lFLOAT  (HOOD) 

GIB  Continue 

C •••••... a. ....a.. •»•••■ a 

c 

C OL  B2I  CREATES  LEPI  ZERO  BUFF  E R 


»C 


IFINS.GT.I19R)  00  TO  200 
40  TO  2vl 

200  PRIM  202 

20?  format  1 • ns  is  tou  l*rof*> 

201  continue 

TLOB  • TINPOIND 
NP01NTS  a NPOINTS-NK"IFT 
TlEFT  a AMN1F IT  111 «T INPOINTSi 1 
IF  ITLEFT.LT. 0.01  A0.50 
PRINT  ISO.  OlSf ANCE.N0I5T0R 
STOP 

TBIG  a AMXJF  (Till  »T  (NPOIMS)  I 
IF  ITBIO.bT.PERIOUl  ftO.TO 
PRINT  1 FiO • OTSTANCf  .NOI5TON 
STOP 

CONTINUF 

IFlAO  • NSN I P T a 1) 

TLOB  a 1111.1,0 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCLCCCtCCCCLLCLCCCCCCCCLCCCCCCCCCCCCCCCL 


AO 


50 


ftO 


TO 


This  LOOP  ADJUSTS  VELOCITIES  AT  1 Mt  SHOCKS  ANO  SHIFTS  OUT  POINTS 
IN  THE  overlap  REOION  that  OCCURS  BREN  SHOCKS  Overtake 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCLCCLCCLCLCCCtCCCCLCCCLCCCCCCCCCCtC 
00  130  INUCXal.NPOTNTS 
NPOlNT  a INOCX-NSMIFT 
TINPOINT)  a T ( I NOE  X I 
UINPOINTI  • KllNOEX) 

IF  ITFLAG.CO.n  GO  TO  |VU 
IF  IT (NPOlNT I .GT .TLOBI  GO  10  120 
TAVG  a (TlNP0INT)«T(HPOlNT-in/2.0 
IFLAfl  a | 

I a NPUINT-2 
GO  TO  90 


1-1 


120 
1)0 
1 • 0 


NSHlFT  a NSHtFT.I 
T ( I •?)  a TINPOINTI 
U I I •? I ■ UINPOINTI 
IF  ITAV8.LT.TIIH  go  TO  80 

u«!*n  b 11 1 1 1 • ii) 1 1 « 1 1 • , 1 1 )•' 1 Avb-i ( n 1/ it ( i*ii-t ( n 1 

Titan  a TAVG 
60  TP  |)0 

IF  ITINPOINTI.LT. TAVOI  GO  10  110 

UINPOINT-ll  • U(NP01NT-n*IU'NP0INTI-UINP01NT-|liPlTAVG-- 
* I NPOTNT* I | )/  IT  INPOINT  1-1  INMOINT-ll » 

TlNPOINT-11  B TAVG 
IFLBfl  a 0 
00  TO  1 20 
NSMIFT  a nShJFT.) 

GO  TO  |)0 

TLOB  • T(NPOTNT) 

CONTINUF 

NPOINTS  a NPOINTS-NSHTFT 
RETURN 


FORMAT  1 1 r* | • •Cxr££pEi)  TIME  RASE  '0  Th|  LEFT*.// . 101. 

••or  stance  ••. tn, s. ioa. is. ix. buena  i ionb*i 
1 format  nm.BEicttoED  time  base  '0  ire  right*,//. io» 

• af)iSTANCEB».t|5.S.10*.TS.)»  ••ITEHAl  IONS*i 
END 


109 


SURROUTINE  TJMEFRE  (u.T.NPOINTS.N,NE»P*.IA,SCAl.NM,ITOT.TNA*,LU,||, 

CSIG.MAF I G .0  T . NBf HC ) 

oi**ension  nui.rUi.mii.iuni 
0Im£N5I0N  hELB  *513> 



c 

C SUB.  TtMFFRE  REnovES  »HC  ZERO  BUFFERS  FROM  AN  INTERMEDIATE 
C BAVEFORM  AND  f-ebFORmS  A FFT.  The  ZERO  BUFFERS  af»E  ALSO 

C CmECRIO  TO  SFl  JF  THE  T STILL  E>1ST. 


• RIOhT  ZERO  BUFFER  NOT  LARGE  ENOUGNA) 


C "L  3AA  »EH0¥ES  The  LEFT  ZERO  BUFFER.  NET  IS  EWUAL  TO  THE 
r Nuhre M OF  ZEROES  REMOVED. 


C OL  SOO  BFNUN0FRS  THE  ELE.  I N Tnf  u AND  T AMRAR1ES  ANU  ha  1 NT  a 1 ns 
C th|  CORRECT  TIME  RELATIONSHIP  BEbTEEN  THE  U ELE. 


C RE  SAMPLE  PROOUCCS  N EOualL t SPACEU  POINTS  FROM  ThE  NPOINTS 
C UNEOUALL*  SPACED  POINTS  OBTAINED  FROM  VAVEPRQP. 


CALL  BESAMPLt  IU.T.N. NPOINTS. OT, MEL** I 

I0P»ISC«G 

Tma»*TIN»-T(H 

DTaTma*/ in-1 ) 


C OL  10  SCALES  THE  FLOATING  PT.  ELE.  OF  U BEFORE  An  INTEGER  EFT 


00  :o  j«i.n 

I*(-)»UI JIASCAL 
10  CONTINUE 


C A FOBARO  INTEGER  FF T IS  TAKEN  UN  IhE  bavEFuRm 


ITOT«IFFTOF(I*,nE*PZ.!OR.ISCfLUI 

£NO 


r 
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LISTINGS  OF  PROGRAM  GRAPH 
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PR06RAM  QRARH 
OlMCNSIO*  MIPI6) 

DIMENSION  fT  1180*1 
DIMENSION  HHOOl 

0 1 ME  NS  I ON  ANP(d,200/»SIOMAfA.2COJ« TAMP(?*0>»*SI6MA<?00> 
01 ME NS ION  XP(55) 

COMMON  VPC ( 2^0  » A) 


PROGRAM  GRAPH  SELECTS  CERTAIN  ►RtOUtNCT  COMPONENTS  FHOM  THE 

SPECTRA  PROM  PROGRAM  C500  AND  PLOTS  The  AMPL  I I UUES  OP  THt 
PREOUENCT  COMPONENTS  ASSOCIATtU  TO  k PARTICULAR  PREUUENCY 
VERSUS  THE  SIGMa  OP  THE  PUMP  OP  tACH  PMfOUENCY  COMPONENT. 

POR  MORE  E<RLANaTION  SEE  AML-TR-77-G.  THE  FOLLOWING  IS  A 
list  or  important  program  parameters. 

AY?  -the  BANOmIDTm  op  the  NUMERICAL  PIlTER. 

SIGN*  -THE  MAXIMUM  SIGMA -2  TO  BE  PLOTTED. 
SISCAL*S1GMA/S  -SCALING  PACTUR  PUR  HON.  AXIS. 

UBGO  -THE  M IN,  VALUE  ON  VERT.  AXIS  IN  OP  RE  .0000?  N/M 
DBENO  -THE  MAX  VALUE  0^  VERT  AXIS  IN  UR  Rt  2«10aa-SN/N 
OHSCAL* I0RPN0-0860) / J-SCALlNG  PACTOR  P OR  VERT.  AXIS. 
LPST  -SPECTRAL  NUMBER  Of  P1HST  PRE.  COMPONENT  PLOT. 
LPENO  -SPECTRAL  NO  OP  LAST  PRC.  COMPONENT  PLOT. 

NO-  THE  NO  OP  SPECTRA  PASSEU  10  GRAPH  By  CSOO. 

MX  - I X. 1 -A ) NO  OP  PRE.  COMP.  FOUND  POM  EACH  SPECTMAL  NO. 


C •••••••••••••••••••••••••••••••••••••••••••••••••••••••• 

READ  (2>  A. PI. B.P2. THETA, DISTANCt, UISINC.NPO.EMAX.no 
PRINT  1000,P1.P?»NO 

1*00  P0RM4TI  • PI  - «.PI0.2t*  Pi  • NO  a *,110) 

P A-P  ?-P  I 
P»«Pf*P 1 
»tr*so. 

SI«MI.M, 

C»ST«i 

|PENO«A 

DBSOmgO. 

OBENOmIOO. 

STSCAL-5I6MX/S. 

DB5CAL ■ 10HFNU-DR8C) /3. 

N(1)«F1 


OL  600S  CALCULATES  THt  CENTER  PRE. 
ThFR  |N  THE  « ARRAY . 


OP  IRE  PTLTtR  ANU  STORES 


DO  600%  JO « SO  • 3 
WtJlm(J/3.ff2 
■ U-1  »*BIJ»-P1 
B«J*1)«MIJ»»P1 

IP  (Mlj.n  .GT.PRaX)  go  TO  6006 
6009  CONTINUE 
6006  LOTaj 


OL  200  HFGINS  THE  FILTERING  OP  EACH  SPECTRUM  MY  READING  IT  PROM 
MAG  TAPr. 


DO  200  ja 1 .NO 
JAVa| 

READ  «2)  SIOpILM.PT 
IP'SIG.GT.SIGMX)  GO  TO  ?01 


DL  100  IS  THE  ACTUAL  NUMERICAL  PlLTEM. 


C«. •••••• •••*•••••••••••••••••••••••••••' 

DO  100  lal.ILM.? 

101  CONTINUE 

CHECK.ABS (B I JAY ) *P  T ( I ) ) 

IPlCMrCA.LT.XY7)  00  TO  203 
IPlPTlII.LT.BIJATl)  GO  TO  100 
JAfajAY* 1 
00  TO  1VI 

203  CHPCA2a*RSlBljAV)-PT M*2) ) 
IPICHECA2.LT.IY7I  GO  TO  2oa 
GO  TO  209 


?U*  IP (LHECA.LT. ChEC«2)  GO  TO  ?06 
GO  TO  luO 

206  IP(PT(I*D.LT.PTII*3))  GO  TO  100 
? 0 a CONTINUt 

If  IJMY.fO.il  GO  TO  21  0 

IF (JAY. |0.?)  GO  TO  211 

IP (JAY. EQ. 3)  00  TO  213 

M*.M4«1 

MRUNmHa 

GO  TO  21A 

210  H l .M 1 « I 
NNltNaMl 
GO  TO  ? 1 A 

211  R?aM?« I 
MRiiNaM2 
GO  TO  214 

213  H 3 aM  3*1 
MRiiN.m3 

214  CONTINUE 

Amp i JAY .MR UN  I aPT  ( I • 1 • 

SIGMA (jAY.MRUN)«S|G 
JAYajAY. | 

If (JM*. to. 91  60  TC  200 
100  CONTINUE 
?00  CONTINUE 
?01  CONTINUE 

c»* •••••••••••••••••••••••••••• •••••♦•••« 





ol  Roi  plots  the  0Est»Eo  spectral  no.  plots. 


C 

DO  801  laLPST.LMENO 
1PM.tO.lt  MffUWvMl 
IP  M. €0.21  MRUN.M2 
lPM.fO.1l  MRUN.MJ 
IP  M. CO. A)  MMUNsMA 
LHa?HOB 
LBL "3hS 1 6 

CALL  PLOT  (6.S.-I2. »«3) 
CALL  ORIGINI I . • 3.0 ' 

CALL  PLTAXIS(0..*.«  S..I 

Call  pltaxtsio.,*..3  .«< 


.......  i a .< 


lO.O.SiGMA.l ..LBL • - ».SI 
i 1 >MGO . UBf NO ,S. .LB .2. A) 


800  TPANSPPPS  The  PARTICULAR  UAIA  OP  ONt  SPECTRAL  NO, 

plotted  to  tme  plotting  arraies. 


MM..M 

TO  BE 


DO  MOO 

VE  T I J « ? I ■ T AMP ( J) "AMP  • I , J> 
VE7 ( J. | 1 .ISIGMA ( J) .SIGMA ( 
800  CONMNUF 


Xh AR. GAMMA. TOL«0. 

MM.MPUN 

IVFC-200 
KKalO 
K3*Kff«l 
A 1 Pa? 

C#. ..... ..a........ ••*#•«•••••••»••*••• •••••••••• ••••••• •* 


PERFORMS  A POLTNOMlNAL  OF  ORDER  **  LEAST  SUdARES  BEST  FIT  ON  THE 
DATA  TO  R£  PLOTTED. 


CALL  PFITLSOIVCC.IVEC.XP.MM.AX.AJ.KIP 
CTSO.TMAX) 

UO  P 1 A jal.MRUN 
YAMP ( J).VEC ( J.9) 

BI A CONTINUE 


XBAR.GAMMA.TOL.TAV.YBAR, 


PRINT  AA4.MRUN 

aaa  Format ( • mrun  . a.iisi 

CALL  PLTOATA  (XSIGMA.VAXPP.MRUN.O. 

call  plot (o.v.o.oo) 

80)  CONTINUE 
End 


0 • U . * S 1 SCAL .UBGO.OB5CAL ,0.08) 
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LISTINGS  OF  PROGRAM  TIMEPIOT 
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I 


1 
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PROORAN  TlRE**L01 
DIMENSION  T («200> .U(*Z00) 

c 

POO€*AR  T|»*fPLOT  REMOVES  The  LEFI  ZERO  BUFFERS  » ROR  SELECTED 
■AVfFORRS  AND  PLOTS  The  FIRST  PURI  ion  OF  THE  SLLtCTEC 
•AVfFOftM S.  POP  NOPE  DC/AJl  SEE  ahl-TN-77-*.  the 
POLL  OP  I NO  IS  A LIST  OF  IrPOMIARI  PROGRAM  PARAMETERS. 

NO-  The  no  OP  PAVEFOPMS  TO  he  LONSIOtatO. 

I HOP*  THE  NO  of  PAPEFORNS  SNIPPED  HilPtEN  Ih£  PLOTTED  ONES  IS 
IMOP-I 

RHOPMmOP  USIO  FOR  INITIALIZATION. 

TNAI* •##?* >46  THE  LCNGTm  OF  The  PLOTTED  PORTION  OF  THE  PAVEFORn 
TSCAL-TRAI/S  SCALING  FACTOR  FUR  HUP.  AMS. 

TTICP.0004F  SPACING  IN  SEC.  ON  HOP.  AMS  TIC  RAPPS. 

4NP.4  MAI.  PARTICLE  VELOCITY  |M  R/S. 

APPSCLP**ARP/?.S  SCALING  FACTOR  FUR  VERT  . AMS. 

IFNOT  USED  to  PRINT  THE  Tic  SPACING  AT  The  0OTTUR  OF  THE  PLOT 
PAPER. 

IFLTT-  USEO  to  POSITION  THE  PLOTS  OR  THE  PLOTTING  PAPER. 


NOMOO 
I HOP  *4 
RHOPhI HOP 
THAI*. OPE*. 44 
TSC»L*TRAA/S. 
TTICp.OPOOZ 


arpscl*E.*arp/2.s 

IfnotpS 

hop-0. 

IFLTT-1 


A ORIGIN  IS  ESTABLISHED  ON  ThE  PLOTTING  PAPER. 

C 

C— 

CALL  PLOTIO..»lf..-J» 

CALL  OPIGlRlO.**.) 

C 

PL  100  IS  The  NAJON  I OOP  OF  THE  PMOORAN.  THE  LUOP  RE AOS  A 

hAVEFORR  FROR  MAG  TAPE  AND  DECIDES  TO  PLOT  IT  UM  TO  READ 
ANOTHER  ONE. 

c 

c •••••• 

00  100  ItPl.NO 

REAO  13)  T.u.N.II.SIO 

IF  (II. EQ. I HOP) 00  TO  999 

90  ro  ioo 

999  CONTINUE 

IhOPpIhOP-hhOP 

C*—«******— 

c 

PL  300  COUNTS  ThE  ZEROES  IN  THE  LEFT  ZtMO  BUFFER.  THE  NO.  OF 
ZEROES  EQUALS  JJ-1. 

C 

C*——***— 

00  300  JJ*).N 

IF  (U  I JJ) .EO.O. ) 00  TO  300 

00  TO  3-0 

io*  continue 

300  CONTINUE 

TSL IPpT ( JJ) 

NRaN-jJ-1 

c 

PL  3)0  RfNOwfS  The  LfFT  ZERO  PUFFER  ANU  SLIPS  The  TINE  AMS 

RACPGAPO  BT  TSlIF  TO  NAINIAIN  The  ORIGINAL  TINE  RELATION- 
SHIP BCThEIn  ThE  ELES.  OF  The  U ARRAY. 

C 

f 

DO  310  iPltNN 
■Ap I • JJ- I 
UltlPUlNA) 

T 1 1 ) aT INA)>TSl IP 
310  CONTINUE 


PL  SSS  SFlGCCTS  The  NO.  OF  u ELES.  NECESSARY  FOR  A PLOT  YNAi 
LONG. 


00  SSS  JPl.N 

IF  fTMAI.LT.Tlj) l 00  TO  400 
SS*  CONTINUE 

400  Continue 

IFflFLTT.EO.l)  RO  TO  • 

IFUFlTT.EO.*)  GO  TO  * 

BO  TO  6 

A CALL  0R10INI T..-T.I 
IFLTT*3 
OD  TO  T 

S CALL  ORIOlNiO.O.S) 

IFLTTP3 
00  TO  T 

4 CALL  0RIBINI4. .).SI 
IFLTT#1 
T CONTINUE 

lBpImahp 

LBLP4HTINE 

CALL  PLTAMM0..*l.ZS.Z.S.«O..*ANP.AMP.|.,LB.3.l  I 
CALL  PLTAlIS(0.*0.*S..O. .0. .TPAl.TTICtLBL.-A.O) 

CALL  PLTO*TAfT.u»J.O«O.O..TSCAi .ANPSCL .0.0B> 

CALL  PLOTLlNEI.3t-l.«..lA*O.I 
PRINT  3J.1I.SI0 

33  FpPRATI*  I NOE ■ * *.IS.*  S I ORA  ■ *.F4.3) 

IF  1 1 PRO T . F 0 . S i no  TO  PR* 

BO  TO  100 

CONTINUE 

iPNOTRlO 

CALL  OlOTLlNt'.3.-|.7..lA.0.> 

PRINT  3* . T T I C 

3*  FOP-AT  IF  10. 4«*  SECONDS  BfTGElN  I JC  RAPPS  ON  HO*.  AMS*I 
IOO  CONTINUE 
IRQ 


APPENDIX  D 


ALTERNATE  REPRESENTATION  OP  THE  IMPROVED  FM  SOLUTION 
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Is 


where 


The  Earnshaw  solution  for  the  two- frequency  boundary  condition, 
u(°,t)  = uop  sin  cot  + Uow  sin  o^t  , (2.10) 

u(x,t')  = uQp  sin  oupT  + u^  sin  o^cp  , (2.11a) 


T = t'  + ^ (uQp  sin  oyp  + uqw  sin  atf>)  . (2.11b) 

Co 

The  weak  signal  term  of  the  Earnshaw  solution  is 

Uw  = uow  sin  “w7  * 0>-l) 

The  improved  FM  solution  (Eq.  2.19,  a time  domain  solution)  was  obtained 
from  Eq.  D.l  by  using  the  following  approximation  of  Eq.  2.11b, 

? = t'  + UQp  sin  co  cp  . (2.16) 

co 


A dimensionless  representation  of  Eqs.  D.l  and  2.16  is  used  to  derive  a 
frequency  domain  representation  for  the  improved  FM  solution.  The  weak 
signal  term  becomes 


where 


V = sin  M 
w 


(D.2) 


V 

w 


ow 


(D.3a) 


and 


$ = y + tfp  sin 


(D.}b) 
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y = <Dpt' 


Equation  D.2  may  also  be  represented  as 


V = Iae1 
w 


if  the  imaginary  part  of  the  result  at  the  end  of  the  derivation  is  taken 
to  equal  V . The  Laplace  transform  is  taken  from  Eq.  D-5, 


r v 

^ Jn  w 


(y)  e"Sy  dy 


The  transform  becomes,  after  one  integration  by  parts  and  a change  from 
an  integration  over  y to  one  over  $, 


V (s)  = 
wv  ' 


-V  (y)  e‘sy 


ri 
J 0 


ft  e"sy  eJM 


With  the  substitution  of  an  expression  for  y obtained  from  Eq.  D.Jb, 
Eq.  D.7  becomes 


rrsT  - i + 3 a T 
w s S*y0 


sa  sin  $ 


A Bessel  function  identity  is  needed  to' evaluate  Eq.  D.8.  A generating 
function  for  the  modified  Bessel  function  of  first  kind,  In,  is 


i z(t+l/t)  « 

. E 

n=-oo 
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10 

Let  t=jed  ; Eq.  D.9  becomes 


z sin  6 


■ E ^ 


(D.10) 


n=-oo 


Equation  D.8  becomes,  with  the  application  of  Eq.  D.10, 


V (s)  = i + V I (a  s)(-j)n  e 

w s s n P Jn 


[-s+(fl+n) ] j<J> 


d«> 


( D . 11 ) 


Evaluation  of  the  integral  in  Eq.  D.ll  is  straightforward  and  yields 


, _2i_  I (a  s) 

tjt)  = i + & y (.j)n  jsLjlL 

w'  ' s s Z_ J v s-j(fi+nj 

n — u ' ' 


(D.12) 


n=-oo 


The  inverse  Laplace  transform  of  Eq.  D.12  is  expressed  as 


i r f-,  °°  _ I (a  s) 

'»w  = SjX  [s  + 30  _E  (-J> 


n=-« 


(fl+n)] 


esy  ds 


(D.13) 


The  integrand  of  Eq.  D.13  has  simple  poles  at  s=0  and  at  s=j(fi+n).  Using 
the  theorem  of  residues  to  evaluate  the  boundary  integral  of  Eq.  D.13 
yields 

V»>-  E . (»•*> 

n=-«> 

where  the  Bessel  function  identity 


In(Jx)  = jn  Jn(x) 


(D.15) 


has  been  used.  Taking  the  imaginary  part  of  Eq.  D.l4,  the  final  result 
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is  obtained 

00 

vw(y)  = 2 FT7  Jn^+n)an  sln(fi+n)y  • (d.i6) 

n=-oo  ^ 

Equation  D.l6  is  an  alternate  representation  for  the  improved  FM  solution 
presented  in  Chapter  II  (Eq.  2.19)  and  is  equivalent  to  Eq.  2.20.  The 
derivation  presented  here  is  due  to  Blackstock.  A similar  approach  was 
used  by  Lockwood^  in  his  derivation  of  a solution  equivalent  to  Fenlon's 
solution. 

Equation  D.l6  shows  an  asymmetry  structure  similar  to  the 
approximate  Fenlon  solution  (Eq.  2.23);  the  argument  of  the  Bessel 
function  is  ft+n,  not  simply  ft  like  in  the  FM  solution.  The  amplitude 
factor  for  Eq.  D.l6  is  ft/(ft+n)  rather  than  an  amplitude  factor  similar  to 
the  one  in  Eq.  2.23-  The  reason  for  the  difference  in  amplitude  factors 
is  that  Eq.  2.23  incorporates  to  some  degree  the  mutual  sideband  formation 
from  both  signals,  while  Eq.  D.l6  incorporates  only  the  sidebands  from 
the  effect  of  the  pump  on  the  weak  signal. 
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